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Integration  of  damage  tolerance  module  with  ASTROS 


FOREWORD 


This  is  the  final  report  on  the  work  performed  by  Knowledge  Systems,  Inc.  on  the  U.S.  Air  Force 
Contract  F33615-96-C-3215,  ”An  ASTROS  Compatible  Strategy  for  Evaluating  the  Aeroelastic 
Response,  Buckling  and  Integrity  of  Composite  A/C”.  This  report  contains  3  parts:  I)  ’’ASTROS 
Damage  Tolerance  Module:  Final  Report”;  2)  ’’ASTROS  Damage  Tolerance  Module:  Interface 
Design  Document”;  and  3)  ’’ASTROS  Damage  Tolerance  Module:  User's  Manual”. 

This  report  details  the  work  performed  to  enhance  the  capability  of  ASTROS  to  perform  prelim¬ 
inary  design  optimization  of  metallic  and  composite  material  aircraft,  based  on  damage  tolerance 
requirements.  The  customized  damage  tolerance  models  that  have  been  implemented  in  ASTROS, 
at  present,  are: 

1 .  Discrete  Source  Damage  Model:  A  lead  crack  in  a  stiffened  panel  with/without  the  presence 
of  a  central  broken  stiffener; 

2.  BuckDel  model:  Buckling  of  a  composite  panel  in  the  presence  of  a  delamination; 

3.  Straight  Crack  Model:  A  panel  with  a  central  crack; 

4.  Rivet  Hole  Crack  Model:  One  (or  two)  crack(s)  emanating  from  one  side  (or  both  sides)  of 
a  rivet  hole; 

5.  Curved  Crack  Model:  A  panel  with  a  curved  crack; 

6.  Rivet  Hole  Curved  Crack  Model:  One  (or  two)  curved  crack(s)  emanating  from  one  side  (or 
both  sides)  of  a  rivet  hole; 

7.  Surface  Crack  Model:  One  centered  surface  crack  in  a  plate; 

8.  Rivet  Hole  Corner  Crack  Model:  Two  comer  cracks  emanating  from  both  the  sides  of  a 
straight-shank  rivet  hole. 

The  authors  acknowledge  the  contributions  of  D.S.  Pipkins,  P.E.  O'  Donoghue,  K.  O'  Sullivan, 
and  H.  Kawai  to  various  parts  of  this  report. 

It  is  a  pleasure  to  acknowledge  the  constant  support,  constructive  criticism,  and  valuable  in¬ 
sights,  provided  by  Drs.  V.A.  Tischler  and  V.B.  Venkayya  of  AFRL  during  the  course  of  this 
project. 
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CHAPTER  I 


INTRODUCTION 


1.1  Introduction 

In  addition  to  requirements  such  as  strength,  stiffness,  and  aeroelastic  response  it  is  necessary  to 
design  both  metallic  and  composite  aircraft  structures  to  withstand  the  effects  of  damage.  By  doing 
so,  safe,  economical  fleets  with  high  operational  readiness  can  be  insured.  The  importance  of  safety 
and  operational  readiness  to  profitability  and  competitiveness  in  a  commercial  aviation  setting 
and  national  defense  in  a  military  one  is  apparent.  To  insure  that  metallic  aircraft  structures  are 
designed  and  maintained  to  withstand  the  effects  of  damage,  the  Federal  Aviation  Administration 
(FAA)  and  United  States  Air  Force  (USAF)  have  established  specific  guidelines  which  must  be 
followed.  The  FAA  requires  commercial  transport  aircraft  certified  under  part  25  of  the  Federal 
Aviation  Regulations  (FAR)  to  meet  certain  Damage  Tolerance  Requirements  (DTR).  Similarly, 
the  USAF  has  specified  a  set  of  DTR  for  metallic  structures  which  are  described  in  detail  in  MIL- 
STD-1530A.  For  composite  structures,  the  FAA  does  not  have  formal  requirements  for  certification 
at  this  time.  The  Air  Force,  on  the  other  hand,  has  specified  a  set  of  DTR  for  composite  structures 
in  AFGS-87221  A.  The  DTR,  as  set  forth  by  both  the  FAA  and  the  Air  Force,  essentially  state  that: 

•  the  residual  strength  of  an  airframe  structural  component  shall  not  drop  below  that  required 
to  sustain  limit  load,  and 

•  that  inspections  must  be  scheduled  to  insure  that  the  required  level  of  residual  strength  is 
maintained. 

The  primary  means  by  which  the  DTR  are  satisfied  by  commercial  and  military  airframe  man¬ 
ufacturers  and  maintenance  organizations  is  through  the  performance  of  a  Damage  Tolerance 
Assessment  (DTA).  The  DTA  involves  the  development  of  damage  growth  curves  and  residual 
strength  diagrams  for  individual  structural  components  of  an  airframe.  This  allows  the  residual 
strength  as  a  function  of  aircraft  usage  to  be  determined.  With  this  information,  manufacturers 
can  design  components  which  minimize  the  evolution  and  growth  of  damage  and  its  effect  on  the 
integrity  of  the  aircraft  structure.  Further,  appropriate  inspection  intervals  can  be  specified  which 
insure  that  the  structural  integrity  of  the  aircraft  will  be  maintained  through  out  its  life. 

At  present,  there  is  no  other  software  which  integrates  damage  tolerance  with  the  other  disci¬ 
plines  (i.e.  strength,  stiffness,  aeroelastic  response,  etc.)  which  impact  the  design  of  an  aircraft 
structure.  As  a  result,  the  design  of  an  aircraft  structure  which  satisfies  damage  tolerance  require¬ 
ments  in  addition  to  those  of  other  disciplines,  is  presently  accomplished  via  a  manual  “cut  and 
try”  procedure.  This  type  of  design  process  is  time  consuming  and  therefore  very  costly  [Nees 
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(1995)].  To  address  this  deficiency  in  existing  software,  this  SBIR  project  is  implementing  a  dam¬ 
age  tolerance  module  into  the  multidisciplinary  analysis  and  design  software,  ASTROS. 

Phase  I  of  this  SBIR  project  has  addressed  the  damage  tolerance  analysis  of  aircraft  structures 
made  of  lanunated  composites.  The  damage  considered  was  in  the  form  of  a  delamination  between 
lamina  in  a  stiffened  panel.  The  Phase  I  effort  accomplished  the  following. 

•  Damage  tolerance  software  applicable  to  laminated  composites,  called  BUCKDEL,  was  de¬ 
veloped.  The  BUCKDEL  software  performs  a  geometrically  nonlinear  analysis  of  stiffened 
laminated  composite  plates  with  and  without  delaminations.  BUCKDEL  allows  the  user  to 
perform:  a  linear  static  solution;  a  linear  buckling  (eigenvalue)  analysis;  and  a  nonlinear 
post-buckling  analysis  through  both  limit  and  bifurcation  points.  BUCKDEL  also  calculates 
the  pointwise  energy  release  rates  around  a  delamination  front  using  an  Equivalent  Domain 
Integral  (EDI)  technique.  The  energy  release  rates  can  be  used  to  predict  the  growth  and 
onset  of  unstable  propagation  of  a  delamination. 

•  The  feasibility  of  using  a  global-local  approach  to  link  an  ASTROS  finite  element  analysis 
(global)  with  a  local  damage  analysis  (such  as  that  performed  by  BUCKDEL)  was  investi¬ 
gated.  This  global-local  approach  was  found  to  be  workable,  due  in  large  part  to  the  AS¬ 
TROS  system  architecture  which  allows  the  user  to  introduce  special  purpose  modules  by 
making  use  of  the  SYSGEN  program  which  is  provided  with  ASTROS. 

Based  on  the  Phase  I  results,  the  implementation  of  a  damage  tolerance  module  into  ASTROS 
which  accounts  for  typical  damage  in  both  composite  and  metallic  structure  is  feasible.  This  report 
describes  the  Phase  II  effort  in  the  implementation  of  the  damage  tolerance  module  in  ASTROS. 
In  addition  to  this  report,  the  users'  manual  and  the  interface  design  document  for  the  damage 
tolerance  module  are  prepared  as  part  of  the  documentation  for  this  project. 

In  phase  II  of  this  SBIR  project,  customized  Damage  Tolerance  Models  (DT  model)  are  im¬ 
plemented  in  ASTROS.  Customized  damage  tolerance  models  are  idealized  damage  scenarios  that 
are  to  be  considered  for  damage  tolerance  requirements.  Localized  damage,  such  as  a  crack  or 
delamination,  is  very  small  in  size  when  compared  to  the  finite  element  model  used  in  preliminary 
design  optimization.  Therefore,  these  damages  are  not  modeled  explicitly  in  the  FEM  model  for 
preliminary  design  optimization.  The  damage  tolerance  module  automatically  generates  detailed 
computational  models  for  damage  tolerance  assessment,  using  a  few  parameters  for  the  definition 
of  the  DT  models  via  bulk  data  cards  introduced  by  the  damage  tolerance  module.  This  feature 
of  the  damage  tolerance  module  greatly  simplifies  the  model  preparation  effort  for  the  user  in  the 
preliminary  design  phase.  The  customized  DT  models  that  have  been  implemented  in  ASTROS  at 
present  are: 

•  Discrete  Source  Damage  Model:  a  lead  crack  in  a  stiffened  panel  with/without  the  presence 
of  a  broken  central  stiffener  [Fig.  2.7] 

•  BuckDel  Model:  buckling  of  a  composite  panel  in  the  presence  of  a  delamination  [Fig.  2.8] 

•  Straight  Crack  Model:  a  panel  with  a  centered  crack  [Fig.  2.1] 


2 


•  Rivet  Hole  Crack  model:  one  (or  two)  crack(s)  emanating  from  one  side  (or  both  sides)  of  a 
rivet  hole  [Fig.  2.2] 

•  Curved  Crack  model:  a  panel  with  a  curved  crack  [Fig.  2.5] 

•  Rivet  Hole  Curved  Crack  model:  one  (or  two)  curved  crack(s)  emanating  from  one  side  (or 
both  sides)  of  a  rivet  hole  [Fig.  2.6] 

•  Surface  crack  model:  one  centered  surface  crack  in  a  plate  [Fig.  2.3] 

•  Rivet  Hole  Corner  Crack  model:  two  corner  cracks  emanating  from  both  side  of  a  straight- 
shank  rivet  hole.  [Fig.  2.4] 

A  master  element  approach  has  been  implemented  in  ASTROS  to  access  these  damage  tol¬ 
erance  models.  Using  these  models  the  user  can  evaluate  fracture  parameters  or  perform  fatigue 
crack  growth  analyses.  Damage  tolerance  based  constraints  for  the  design  optimization  can  then 
be  used  in  the  preliminary  design  optimization. 

1.2  SBIR  Technical  Objectives 

The  objective  of  this  SBIR  project  is  to  develop  damage  tolerance  software  as  an  analysis  mod¬ 
ule  in  the  multidisciplinary  analysis  and  design  software,  ASTROS.  The  software  uses  state  of 
the  art,  computationally  efficient  algorithms  for  determining  the  residual  strength  and  life  [Atluri 
(1995)]  of  metallic  and  composite  structures  with  damage.  It  enhances  the  existing  capabilities 
of  ASTROS  by  allowing  constraints  based  on  damage  tolerance  requirements  to  be  considered 
simultaneously  with  those  based  on  strength,  stiffness,  aeroelastic  response,  etc.  during  design 
optimization.  Included  in  potential  commercial  applications  of  such  a  capability  are  industries 
involved  in  the  design  of  aircraft  structures,  automobiles,  bridges  and  buildings. 

ASTROS  is  well  suited  for  modeling  the  global  strength,  stiffness,  and  aeroelastic  response  of 
an  undamaged,  stiffened  structure.  However,  it  presently  does  not  have  the  capability  to  account  for 
local  damage  such  as  cracks,  delaminations,  and  penetration  holes.  Therefore,  this  SBIR  project  is 
to  develop  a  damage  tolerance  module  for  ASTROS.  The  damage  tolerance  module 

•  utilizes  the  existing  capabilities  of  ASTROS  for  performing  the  modeling  of  the  structure; 

•  generates  load  spectra  in  terms  of  ASTROS  load  cases; 

•  develops  customized  damage  tolerance  modules  to  model  local  damages  in  metallic  and 
composite  materials; 

•  has  a  seamless  interface  to  the  other  modules  in  ASTROS,  via  an  engineering  databa.se  and 
high  level  programming  language  (MAPOL). 

•  defines  new  bulk  data  entries,  relational  schema  and  error  messages  needed  to  provide  a 
consistent  interface  as  the  other  modules  in  ASTROS. 

The  damage  scenarios  considered  in  the  damage  tolerance  module  are: 
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Figure  1.1:  Central  Crack  and  Broken  Stiffener  in  a  Panel 


•  single  or  multiple  skin  cracks  (Widespread  Fatigue  Damage)  in  a  stiffened  structure,  in¬ 
cluding  the  effect  of  broken  stiffeners  (Fig.  1.1)  [tensile  loading,  metallic  and  composite 
materials]; 

•  skin  crack  turning  at  stiffeners  [tensile  loading,  metallic  and  composite  materials]; 

•  single  or  multiple  part  ellipticaFcircular  surface  flaws  at  holes  or  other  stress  raisers  (Fig.  1 .2) 
[tensile  loading;  metallic  materials] 

•  holes  (Fig.  1.3)  [compressive  loading,  metallic  and  composite  materials];  and 

•  delaminations  (Fig.  1.4) [compressive  loading,  composite  materials]. 

For  metallic  or  composite  aircraft  structures  loaded  in  tension,  the  damage  of  principal  interest 
during  design  is  usually  in  the  form  of  cracks  normal  to  the  direction  of  principal  tension.  These 
cracks  typically  occur  over  time  due  to  fatigue  or  suddenly  due  to  an  event  such  as  an  uncon¬ 
tained  failure  of  a  rotating  engine  component  or  battle  damage.  A  crack  arising  suddenly  due  to  a 
catastrophic  event  is  an  example  of  discrete  source  damage  (DSD).  In  reality,  DSD  in  a  structure 
loaded  in  tension  such  as  a  fuselage  or  lower  wing  would  be  in  the  form  of  an  irregular  shaped 
hole.  However,  since  a  crack  of  length  'a'  is  more  critical  than  a  hole  of  diameter  'a' ,  the  crack 
representation  of  DSD  represents  a  worst  case  scenario.  Therefore,  for  reasons  of  conservatism 
and  practicality  (both  experimental  and  analytical)  the  crack  is  used  in  the  certification  of  primary 
structural  elements  loaded  in  tension.  Residual  strength  and  life  calculations  for  structures  con¬ 
taining  cracks  and  loaded  in  tension  will  be  based  on  Linear  Elastic  Fracture  Mechanics  (LEFM). 
Thus,  the  parts  of  the  damage  tolerance  modules  which  analyze  cracked  structure  will  compute 
stress  intensity  factors  and  their  sensitivities  to  changes  in  the  design  variables.  These  values  can 
then  be  used  to  evaluate  constraints  based  on  residual  strength  requirements.  For  constraints  based 
on  residual  life  (fatigue)  requirements,  the  computed  stress  intensity  factors  will  be  used  in  crack 
growth  equations  to  determine  the  time  required  for  a  crack  or  cracks  to  grow  to  a  critical  size.  The 
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Figure  1.2:  Failed  Lower  Wing  Panel  of  a  U.S.  Air  Force  C-141B  Due  to  the  Growth  of  a  Surface 
Flaw 


Figure  1.3:  Hole  in  the  Upper  (Compression)  Skin  of  a  Wing 


Figure  1.4:  Delamination  Due  to  Impact  of  a  Laminated  Composite 


critical  size  is  determined  from  a  residual  strength  requirement  such  as  that  requiring  the  structure 
to  be  able  to  sustain  limit  load  at  any  time  in  its  life. 

For  metallic  or  composite  aircraft  structures  loaded  in  compression,  the  primary  concern  is  with 
DSD  in  the  form  of  a  hole  or  crack  parallel  to  the  direction  of  principal  compression.  For  laminated 
composites,  delaminations  between  lamina  are  also  of  concern.  Delaminations  can  result  from 
excessive  interlaminar  shear  stresses  or  through-the-thickness  tensile  stresses  at  holes,  free  edges, 
section  changes,  or  in  bonded  joints;  panel  buckling;  and  low  velocity  impact.  The  modeling  of 
delaminations  in  laminated  composites  was  addressed  in  Phase  I  and  resulted  in  the  development  of 
the  BUCKDEL  software.  The  residual  strength  of  structures  loaded  in  compression  will  be  based 
on  the  buckling  and  post-buckling  behavior  of  the  structure.  For  laminated  composites  containing 
delaminations,  the  pointwise  energy  release  rate  around  the  delamination  front  will  also  be  used  in 
the  residual  strength  prediction. 

The  motivation  for  modeling  DSD  in  a  structure  loaded  in  compression  in  the  form  of  a  hole 
or  crack  parallel  to  the  direction  of  principal  compression  is  as  follows.  For  a  stiffened  structure, 
the  buckling  load  is  expected  to  vary  significantly  with  the  size,  shape  (i.e.  circular  or  elliptical), 
and  location  (i.e.  distance  from  the  wing  tip)  of  the  hole.  Some  results  reported  in  the  literature 
[Vellaichamy  et.  al  (1990),  Nemeth  (1990),  and  Britt  (1994)]  indicate,  as  is  to  be  expected,  that  for 
an  elliptical  hole  with  the  major  axis  along  the  direction  of  compression,  the  initial  buckling  load 
is  lower  than  that  for  a  circular  hole  of  the  same  area.  Similarly,  if  a  panel  is  subject  to  pure  shear 
in  the  x-y  plane,  the  shear  buckling  load  will  be  minimum  when  the  major  axis  of  the  ellipse  is  at 
45  degrees  to  the  x  axis.  The  buckling  load  will  decrease  with  an  increase  in  the  aspect  ratio  of  the 
ellipse.  The  most  severe  case  being  in  the  limit  when  the  ellipse  degenerates  into  a  crack  (with  the 
crack  axis  along  the  direction  parallel  to  the  direction  of  compression).  Based  upon  this  discussion, 
it  is  anticipated  that  the  worst  case  DSD  scenario  for  primary  structure  loaded  in  compression  will 
be  a  crack  oriented  such  that  its  axis  is  parallel  to  the  direction  of  maximum  compressive  stress. 
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CHAPTER  II 


TECHNICAL  DESCRIPTION 


This  chapter  describes  the  technical  details  in  the  damage  tolerance  module  in  ASTROS.  The 
damage  tolerance  module  in  ASTROS  uses  a  master-element  approach  to  extract  loading  condi¬ 
tions  for  the  damage  tolerance  models.  An  automated  global-local  analyzer  is  used  to  perform 
hierarchical  analysis  of  discrete  source  damage.  An  integrated  geometry  modeler  and  mesh  gen¬ 
erator  is  used  to  construct  damage  tolerance  models  using  a  few  user  specified  parameters.  Finite 
element  alternating  methods  are  used  for  the  evaluations  of  stress  intensity  factors  for  straight  or 
curved  cracks,  3D  cracks  of  elliptical/circular  shape.  BUCKDEL  is  used  for  the  evaluation  of 
the  buckling  load  of  a  composite  panel  with  delamination.  A  load  spectrum  module  is  used  to 
generate  fatigue  loading  history  in  terms  of  ASTROS  loading  cases.  Details  of  these  techniques 
are  described  in  this  chapter.  The  reader  is  assumed  to  be  familiar  with  ASTROS  (Automated 
STRuctural  Optimization  System). 


2.1  Master-element  approach 

A  master-element  approach  is  used  to  enable  the  automatic  generation  of  loading  conditions  for 
damage  tolerance  models.  Localized  damages,  such  as  a  crack  or  delamination,  are  very  small  in 
size  when  compared  to  the  overall  finite  element  model  used  in  the  preliminary  design  optimiza¬ 
tion.  Therefore,  these  damages  are  not  modeled  explicitly  using  FEM  mesh  during  the  preliminary 
design  optimization.  For  each  of  the  damages  specified  in  the  preliminary  design  model,  a  separate 
damage  tolerance  model  is  generated  for  damage  tolerance  assessment. 

A  master  element  is  the  element  in  the  Preliminary  Design  model  (PD  model)  that  contains 
the  damage  under  consideration.  Since  it  is  small  in  size,  the  damage's  effect  on  the  load  flow 
in  the  PD  model  is  practically  negligible.  Therefore,  the  loading  condition  in  the  master  element 
of  the  PD  model  (without  damage)  represents  the  loading  condition  for  the  panel  with  a  damage. 
Based  on  this  assumption  the  damage  tolerance  module  extracts  the  loading  information  from  the 
master  element  from  the  engineering  database  after  the  static  mechanical  analysis,  and  generates 
the  Damage  Tolerance  model  (DT  model)  for  damage  tolerance  assessment. 

The  damage  tolerance  module  automatically  generates  detailed  computational  models  for  dam¬ 
age  tolerance  assessment,  using  a  few  parameters  for  the  definition  of  the  DT  models  via  bulk  data 
cards  introduced  by  the  damage  tolerance  module.  This  feature  of  the  damage  tolerance  module 
greatly  simplifies  the  model  preparation  effort  for  the  user  in  the  preliminary  design  phase. 

There  are  a  number  of  advantages  in  using  the  master-element  approach.  The  master-element 
approach,  together  with  the  high  level  of  abstraction  of  the  damage  tolerance  model,  makes  it 
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possible  for  the  development  of  interchangeable  modules  for  performing  the  damage  tolerance 
assessment  for  the  same  damage  using  different  methods. 

The  master-element  approach  allows  the  damage  tolerance  module  to  treat  different  types  of 
damage  in  a  similar  fashion. 

When  the  master-element  approach  is  used,  the  damage  tolerance  analysis  imposes  no  addi¬ 
tional  requirements  on  the  preparation  of  the  PD  model.  A  single  PD  model  can  be  used  to  assess 
multiple  damage  scenarios  with  a  single  or  multiple  damage. 

However,  the  master-element  approach  can  not  be  used  when  the  size  of  a  damage  is  so  large 
that  the  loading  conditions  around  the  damage  are  altered  significantly  due  to  its  presence.  This 
may  occur  in  the  analyses  of  the  structural  details  where  the  size  of  the  damage  is  not  small  when 
compared  to  the  structural  details  in  the  FEM  under  consideration.  In  such  a  case  the  effect  of  the 
damage  on  the  load  distribution  must  be  taken  into  account  in  order  to  assess  the  capacity  of  the 
damage  tolerance  of  the  structure. 

2.2  Damage  tolerance  models 

Although  the  structural  details  in  the  vicinity  of  a  crack  have  significant  influence  on  the  stress 
intensity  factors  of  the  crack,  such  details  are  to  be  determined  in  the  later  stage  of  design.  There¬ 
fore,  accurate  modeling  of  structural  details  in  the  vicinity  of  a  crack  for  specific  location  is  not 
necessary  in  the  phase  of  preliminary  design.  However,  it  is  necessary  to  model  the  “average”  struc¬ 
tural  behaviors  so  that  the  overall  damage  tolerance  capability  of  the  structure  can  be  achieved  in 
preliminary  design,  where  the  interaction  between  multiple  disciplines  requires  intensive  computa¬ 
tion.  Once  the  preliminary  design  ensures  that  the  structure  has  adequate  overall  damage  tolerance 
capability,  the  detailed  design  (which  does  not  involves  multiple  disciplines)  can  easily  meet  the 
damage  tolerance  requirement  by  changing  the  structural  details  near  the  critical  locations. 

This  section  describes  the  customized  damage  tolerance  models  that  have  been  implemented 
into  ASTROS. 

2.2.1  2D  through  wall  crack 

In  this  damage  model  [see  Fig.  2.1],  a  small  crack  is  located  at  the  center  of  a  panel.  When  the 
dimension  of  the  panel  is  not  explicitly  specified  by  the  user,  the  damage  tolerance  module  will 
calculate  the  dimension  based  on  the  size  of  the  master  element.  A  square  panel  will  be  generated 
with  an  area  equal  to  that  of  the  master  element.  A  coordinate  transformation  is  performed  by  the 
damage  tolerance  module  so  that  the  crack  is  on  the  x-axis  in  the  damage  tolerance  model  as  seen 
in  Fig.  2.1.  The  loading  condition  is  taken  from  the  corresponding  master  element.  A  2D  plane 
stress  analysis  is  performed.  The  beta  factors  for  mode  I  and  mode  II  SIFs  are  defined  as: 

_  Kj] 

\/ tta  ^j^yy^TZa 

shear  stresses  in  the  master  element  in  the  crack  coordinate 

system. 


Pi 


where  Ov  and  Gxy  are  the  normal  and 
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Figure  2. 1 :  A  through  wall  crack 


The  user  can  chose  from  one  of  the  following  methods  to  perform  the  fracture  analysis.  They 
are  i)  2D  finite  element  alternating  method  (FEAM);  ii)  using  user  supplied  beta  factors;  iii)  using 
infinite  body  solution.  In  practice,  one  can  use  the  2D  FEAM  to  calculate  the  beta  factor;  then,  use 
the  beta  factor  in  the  subsequent  analyses. 

2.2.2  2D  through  wall  cracks  emanating  from  a  hole 

In  this  damage  model  [see  Fig.  2.2],  one  (or  two)  small  crack(s)  emanates  from  a  rivet  hole  at 
the  center  of  a  panel.  When  the  dimension  of  the  panel  is  not  explicitly  specified  by  the  user,  the 
damage  tolerance  module  will  calculate  the  dimension  based  on  the  size  of  the  master  element.  A 
square  panel  will  be  generated  with  an  area  equal  to  that  of  the  master  element.  A  coordinate  trans¬ 
formation  is  performed  by  the  damage  tolerance  module  so  that  the  cracks  are  on  the  x-axis  in  the 
damage  tolerance  model  as  seen  in  Fig.  2.2.  The  loading  condition  is  taken  from  the  corresponding 
master  element.  A  2D  plane  stress  analysis  is  performed.  The  beta  factors  for  mode  I  and  mode  II 
SIFs  are  defined  as: 

R  -—^1—  R 

”  (5y^  ’ 

where  <5y  and  (5xy  are  the  normal  and  shear  stresses  in  the  master  element  in  the  crack  coordinate 
system;  a  =  r-f  (ui  +02)12- 


Figure  2.2:  Two  cracks  emanating  from  a  hole 
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Figure  2.3:  A  surface  flaw  of  the  shape  of  semi-ellipse/circle 


When  there  are  two  cracks  emanating  from  the  hole,  these  two  cracks  must  be  aligned  in  a  line 
passing  through  the  center  of  the  hole. 

The  user  can  choose  from  one  of  the  following  methods  to  perform  the  fracture  analysis.  They 
are  i)  2D  finite  element  alternating  method  (FEAM);  ii)  using  user  supplied  beta  factors;  iii)  using 
infinite  body  solution.  In  practice,  one  can  use  the  2D  FEAM  to  calculate  the  beta  factor;  then,  use 
the  beta  factor  in  the  subsequent  analyses. 


2.2.3  Surface  flaw 

In  this  damage  model  [see  Fig.  2.3],  a  surface  flaw  is  located  at  the  center  of  a  panel.  The  surface 
flaw  is  of  a  half  elliptical/circular  shape.  The  major  (or  the  minor)  axis  aligns  with  the  surface; 
while  the  other  axis  is  in  the  thickness  direction.  When  the  length  and  width  of  the  panel  are  not 
explicitly  specified  by  the  user,  the  damage  tolerance  module  will  calculate  the  dimensions  based 
on  the  size  of  the  master  element.  A  square  panel  will  be  generated  with  area  equal  to  that  of  the 
master  element.  If  the  user  does  not  specify  the  thickness  explicitly,  the  thickness  of  the  master 
element  is  used.  A  coordinate  transformation  is  performed  by  the  damage  tolerance  module  so  that 
the  surface  flaw  is  in  the  x-z  plane  in  the  damage  tolerance  model  as  seen  in  Fig.  2.3.  The  loading 
condition  is  taken  from  the  corresponding  master  element.  The  beta  factors  for  mode  I,  mode  II, 
and  mode  III  SIFs  are  defined  as: 


Pi  = 


h  = 


Kii 


Ps  = 


Kni 


where  <3y  and  <5xy  are  the  normal  and  shear  stresses  in  the  master  element  in  the  crack  coordinate 
system;  a  =  min(a,^?). 

The  user  can  choose  from  one  of  the  following  methods  to  perform  the  fracture  analysis.  They 
are  i)  3D  finite  element  alternating  method  (FEAM);  ii)  using  user  supplied  beta  factors.  In  prac¬ 
tice,  one  can  use  the  3D  FEAM  to  calculate  the  beta  factor;  then,  use  the  beta  factor  in  the  subse¬ 
quent  analyses. 


10 


Figure  2.4:  Two  corner  cracks  emanating  from  a  hole 


2.2.4  Corner  cracks  emanating  from  a  hole 

In  this  damage  model  [see  Fig.  2.4],  two  comer  cracks  symmetrically  emanate  from  a  hole  at  the 
center  of  a  panel.  The  two  corner  cracks  have  the  same  size.  They  are  of  a  quarter  elliptical/circular 
shape.  The  major  (or  the  minor)  axis  aligns  with  the  surface;  while  the  other  axis  is  in  the  thickness 
direction  along  the  surface  of  the  hole.  When  the  length  and  width  of  the  panel  are  not  explicitly 
specified  by  the  user,  the  damage  tolerance  module  will  calculate  the  dimension  based  on  the  size 
of  the  master  element.  A  square  panel  will  be  generated  with  an  area  equal  to  that  of  the  master 
element.  If  the  user  does  not  specify  the  thickness  explicitly,  the  thickness  of  the  master  element 
is  used.  A  coordinate  transformation  is  performed  by  the  damage  tolerance  module  so  that  the 
surface  flaw  is  in  the  x-z  plane  in  the  damage  tolerance  model  as  seen  in  Fig.  2.4.  The  loading 
condition  is  taken  from  the  corresponding  master  element.  The  beta  factors  for  mode  I,  mode  II, 
and  mode  III  SIFs  are  defined  as: 


Pi  = 


Ki 


P3  = 


Km 


where  Oy  and  Gxy  are  the  normal  and  shear  stresses  in  the  master  element  in  the  crack  coordinate 
system;  a  =  min(a,^?). 

The  user  can  choose  from  one  of  the  following  methods  to  perform  the  fracture  analysis.  They 
are  i)  3D  finite  element  alternating  method  (FEAM);  ii)  using  user  supplied  beta  factors.  In  prac¬ 
tice,  one  can  use  the  3D  FEAM  to  calculate  the  beta  factor;  then,  use  the  beta  factor  in  the  subse¬ 
quent  analyses. 


2.2.5  2D  through  wall  curved  crack 

In  this  damage  model  [see  Fig.  2.5],  a  curved  crack  is  located  in  a  panel.  The  crack  is  defined 
piecewise-linearly  by  a  number  of  points  on  the  crack.  When  the  dimension  of  the  panel  is  not 
explicitly  specified  by  the  user,  the  damage  tolerance  module  will  calculate  the  dimension  based 
on  the  size  of  the  master  element.  A  square  panel  will  be  generated  with  an  area  equal  to  that  of 
the  master  element.  The  loading  condition  is  taken  from  the  corresponding  master  element.  A  2D 
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Figure  2.5:  A  curved  crack 


plane  stress  analysis  is  performed.  The  beta  factors  for  mode  I  and  mode  II  SIFs  are  defined  as: 

ft  ft  - 

V  ^xy  y 

where  a^,  and  Oxy  are  the  normal  and  shear  stresses  in  the  master  element  in  the  coordinate  system 
for  the  equivalent  crack.  The  equivalent  crack  is  defined  by  the  straight  line  connecting  the  two 
crack  tips  of  the  curved  crack,  a  is  the  half  length  of  the  equivalent  crack. 

The  user  can  choose  from  one  of  the  following  methods  to  perform  the  fracture  analysis.  They 
are  i)  2D  distributed-dislocation-based  finite  element  alternating  method  (DFEAM);  ii)  using  user 
supplied  beta  factors.  In  practice,  one  can  use  the  2D  DFEAM  to  calculate  the  beta  factor;  then, 
use  the  beta  factor  in  the  subsequent  analyses. 

2.2.6  2D  through  wall  curved  cracks  emanating  from  a  hole 

In  this  damage  model  [see  Fig.  2.6],  one  (or  two)  curved  crack(s)  emanates  from  a  hole  at  the 
center  of  a  panel.  The  crack  is  defined  piecewise-linearly  by  a  number  of  points  on  the  crack. 
When  the  dimension  of  the  panel  is  not  explicitly  specified  by  the  user,  the  damage  tolerance 
module  will  calculate  the  dimension  based  on  the  size  of  the  master  element.  A  square  panel  will 
be  generated  with  an  area  equal  to  that  of  the  master  element.  The  loading  condition  is  taken  from 
the  corresponding  master  element.  A  2D  plane  stress  analysis  is  performed.  The  beta  factors  for 
mode  I  and  mode  II  SIFs  are  defined  as: 

a  _  p  Kii 

HI  -  ^  AT-  >  P2  =  - - 7= 

C5y  yj  Ttd  ^xy  \ 

where  <5y  and  Oxy  are  the  normal  and  shear  stresses  in  the  master  element  in  the  coordinate  system 
for  the  equivalent  crack.  The  equivalent  crack  is  defined  by  the  straight  line  connecting  the  two 
crack  tips  of  the  curved  crack,  a  is  the  half  length  of  the  equivalent  crack. 

The  user  can  choose  from  one  of  the  following  methods  to  perform  the  fracture  analysis.  They 
are  i)  2D  distributed-dislocation-based  finite  element  alternating  method  (DFEAM);  ii)  using  user 
supplied  beta  factors.  In  practice,  one  can  use  the  2D  DFEAM  to  calculate  the  beta  factor;  then, 
use  the  beta  factor  in  the  subsequent  analyses. 
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Figure  2.6:  The  curved  cracks  emanating  from  a  rivet  hole 


Figure  2.7:  A  lead  crack  in  a  stiffened  panel  with/without  a  broken  central  stiffener 


2.2.7  Discrete  source  damage 

In  this  damage  model  [see  Fig.  2.7],  one  crack  is  located  at  the  center  of  a  stiffened  panel.  It  can 
be  a  pressurized  curved  panel  or  a  flat  panel.  The  crack  is  either  parallel  to  the  frame  or  to  the 
stringer.  The  user  can  specify  whether  the  center  stiffener,  passing  across  the  crack,  is  broken  or 
not.  The  user  must  specify  the  panel  size  so  that  a  global  damage  tolerance  model  can  be  analyzed 
to  capture  the  load  flow  in  the  panel  due  to  the  crack.  The  loading  condition  for  the  global  damage 
tolerance  model  is  obtained  from  the  master  element.  The  automatic  global-local  analyzer  extracts 
a  local  model  with  boundary  conditions  obtained  from  the  analysis  of  the  global  damage  tolerance 
model.  A  local  analysis  based  on  FEAM  is  performed  to  obtain  the  stress  intensity  factors.  The 
beta  factors  for  mode  I  and  mode  II  SIFs  are  defined  as; 

R  -  R  = 

“  Cy,/^  ’ 

where  Cy  and  Gxy  are  the  normal  and  shear  stresses  in  the  master  element  in  the  crack  coordinate 
system. 
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Figure  2.8:  Buckling  of  a  composite  panel  with  elliptical/circular  delamination 


2.2.8  BUCKDEL 

In  this  damage  model  [see  Fig.  2.8],  an  elliptical/circular  delamination  is  located  at  the  center  of  a 
panel.  When  the  dimension  of  the  panel  is  not  explicitly  specified  by  the  user,  the  damage  tolerance 
module  will  calculate  the  dimensions  based  on  the  size  of  the  master  element.  A  square  panel  will 
be  generated  with  an  area  equal  to  that  of  the  master  element.  A  coordinate  transformation  is 
performed  by  the  damage  tolerance  module  so  that  the  major  axis  of  the  delamination  is  in  the 
jc-axis  direction  as  seen  in  Fig.  2.8.  The  loading  condition  is  taken  from  the  corresponding  master 
element.  The  buckling  load  of  the  panel  is  computed  using  the  BUCKDEL  module. 

2.3  Automated  global-local  analyzer 

The  damage  tolerance  analysis  of  discrete  source  damage  in  aging  aircraft  poses  many  significant 
technological  challenges.  As  fracture  mechanics  treats  essentially  local  phenomena,  very  detailed 
modeling  of  aircraft  structures  is  necessary  to  analyze  discrete  source  damage.  Not  only  local 
features  (such  as  fasteners),  but  also  global  features  (such  as  stiffeners  and  joint  configurations) 
must  be  considered.  The  fact  that  an  aircraft  structure  at  the  global  level  is  a  built-up  structure 
assembled  from  a  large  number  of  parts  complicates  the  analysis. 

The  traditional  approach  to  this  problem  is  to  utilize  a  very  large  number  of  elements  for  em¬ 
bedded  local  details  with  aggressive  mesh  refinement  in  a  coarser  global  model.  The  problem  is 
then  solved  in  a  single  analysis.  This  approach  takes  many  hours  of  a  very  powerful  computer 
to  calculate  reasonably  accurate  results  using  a  modern  workstation  even  in  the  case  of  a  single 
linear  elastic  analysis.  A  parametric  survey  or  design  optimization  can  only  be  achieved  with  a 
supercomputer.  In  addition,  it  is  very  tedious  and  time  consuming  to  create  a  finite  element  mesh 
with  these  global  and  local  features  in  a  single  model.  This  makes  parametric  study  and  design 
optimization,  which  require  drastic  modification  of  the  mesh,  virtually  impossible. 

Generally,  a  multi-stage  hierarchical  approach  is  a  more  efficient  method  for  such  analyses. 
However,  the  inherent  complexity  of  hierarchical  analysis  and  the  absence  of  conventional  pre- 
and  post-  processors  which  directly  support  hierarchical  analyses,  made  this  approach  impractical 
until  recently.  There  are  several  inherent  complexities  in  hierarchical  analysis.  The  first  is  the 
simplification  of  the  structure  at  each  coarse  stage.  The  second  is  the  extraction  of  a  subregion 


14 


for  each  detailed  analysis.  The  third  is  the  transfer  of  boundary  conditions,  i.e.,  the  conversion 
from  the  analysis  result  in  the  previous  coarse  stage  to  the  boundary  of  the  subregion  in  the  next 
detailed  stage.  It  is  obvious  that  sophisticated  support  of  pre-  and  post-processing  is  necessary  to 
circumvent  these  difficulties. 

The  automated  global-local  analyzer  is  based  on  feature  modeling  technology,  which  has  been 
recently  adopted  in  the  field  of  CAD.  Rather  than  dealing  with  a  finite  element  model  directly  at 
each  stage,  it  deals  with  a  feature  model  -  a  high  level  geometry  model.  After  the  feature  model  and 
associated  analysis  conditions  are  defined,  an  arbitrary  number  of  stages  of  hierarchical  analyses 
(including  simplification,  subregion  extraction  and  boundary  condition  transfer)  can  be  performed 
fully  automatically.  Because  the  feature  model  is  parametric  and  defined  by  several  key  design 
parameters,  a  parametric  study  can  be  performed  by  simply  changing  the  values  of  these  design 
parameters. 

2.3.1  Hierarchical  analysis  methodology 

An  automated  global-local  analysis  methodology  is  used  for  the  analysis  of  a  stiffened  structure 
with  discrete  source  damage.*  Since  the  size  of  the  discrete  source  damage  is  usually  in  the  same 
scale  as  a  stiffened  panel,  it  is  necessary  to  model  the  stiffeners  as  well  as  the  skin  in  the  damage 
tolerance  model.  To  reduce  computational  cost  a  global-local  analysis  methodology  is  used.  In 
the  global  model  cracks  are  modeled  via  unconnected  nodes  at  the  crack  locations.  This  crude 
representation  of  the  crack  in  the  global  model  reflects  the  loss  of  stiffness  in  the  structure,  so  that 
the  redistribution  of  loads  among  the  skin,  stringers  and  ribs  can  be  captured.  Broken  stringers 
and  ribs,  if  any,  are  also  accounted  for  in  a  similar  manner.  The  details  of  the  crack  tip  fields 
are  ignored  at  this  level  of  analysis.  The  global  analysis  results  are  used  to  construct  a  free-body 
diagram  (Fig.  2.1 1)  of  the  cracked  sheet  (local  model),  with  the  applied  loading  on  the  sheet  being 
the  reaction  forces  from  stringers  and  ribs  as  well  as  the  loading  on  the  periphery  of  the  sheet.  The 
local  model  is  analyzed  by  the  damage  tolerance  module  using  a  finite  element  alternating  method; 
while  the  global  model  is  analyzed  by  the  existing  ASTROS  module  for  static  analysis. 

To  illustrate  this  methodology,  consider  the  case  of  a  wing  containing  a  crack  in  the  skin  of 
the  lower  surface  [Fig.  2.9].  (Note  that  the  methodology  is  general  in  nature  and  can  be  used  for 
a  structure  containing  damage  in  the  form  of  penetration  holes  and/or  delaminations.)  A  coarse 
finite  element  mesh  is  first  used  to  model  the  global  behavior  of  the  cracked  panel  [Fig.  2.10].  The 
boundary  condition  for  the  global  DT  model  is  obtained  from  the  master  element  in  the  design 
model  using  the  master  element  approach  discussed  in  §  2.1.  The  traction  and/or  displacement 
boundary  conditions  to  be  applied  to  the  local  model  are  determined  from  the  results  of  the  global 
analysis.  These  boundary  conditions  include  the  reaction  forces,  exerted  by  the  stringers  and  ribs, 
on  the  skin.  The  Finite  Element  Alternating  Method  (FEAM)  is  used  in  the  local  analysis  to 
determine  the  stress  intensity  factors.  The  FEAM  allows  a  coarse  finite  element  mesh  to  be  used 

*  Here  the  global  and  local  model  refer  to  the  damage  tolerance  model  generated  by  the  damage  tolerance  module, 

as  oppose  to  the  preliminary  design  model  used  in  the  process  of  design  optimization. 
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Upper  surface 


Figure  2.9:  A  design  optimization  model 


for  the  local  model  of  the  skin  because  the  crack  tip  fields  are  captured  by  an  analytical  solution 
and  thus  cracks  need  not  be  modeled  explicitly  in  the  mesh.  The  FEAM  along  with  the  other 
local  damage  modeling  methodologies  to  be  implemented  in  the  damage  tolerance  module  will  be 
discussed  in  detail  in  later  sections. 

Damage  in  the  form  of  delaminations  typically  need  not  be  accounted  for  in  the  global  model. 
On  the  other  hand,  damage  in  the  form  of  cracks  or  penetration  holes  does  need  to  be  accounted 
for  if  it  is  of  sufficient  size  to  affect  the  load  transfer  in  the  structure. 

In  the  global  analysis,  the  stringers  are  modeled  as  beams  (ASTROS  GEAR  elements)  attached 
to  the  skin  (ASTROS  CQUAD4),  and  ribs  are  modeled  as  plates  (ASTROS  CQUAD4).  To  account 
for  fastener  flexibility,  beam  elements  are  used  to  model  rivet  connections  between  stiffening  el¬ 
ements  and  the  skin.  The  global  model  is  automatically  generated  by  the  integrated  geometry 
modeler  and  mesh  generator  in  the  damage  tolerance  module.  Most  details  about  the  automated 
global-local  analyzer  are  presented  in  the  rest  of  this  section. 

2.3.2  Geometry  Modeling  and  Mesh  Generation 

This  section  describes  the  most  important  ingredient  of  the  automatic  global-local  analyzer:  the 
integrated  GEOmetry  modeler  and  MESH  generator  (GEOMESH).  GEOMESH  generates  the 
global  damage  tolerance  model,  using  the  user  specified  geometry  parameters  for  the  discrete 
source  damage  and  the  loading  condition  in  the  corresponding  master  element  from  the  prelim¬ 
inary  design  model.  It  invokes  a  separate  ASTROS  process  to  perform  the  global  analysis.  After 
the  global  analysis,  it  1)  imports  the  global  analysis  results;  2)  determines  the  area  for  local  analy¬ 
sis;  3)  extracts  a  detailed  boundary  condition  for  the  local  model  (see  the  free  body  loading  diagram 
in  Fig.  2. 1 1).  A  summary  of  how  GEOMESH  works  follows. 
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Figure  2.10:  A  Global  Analysis 


Figure  2.1 1:  The  Isolated  Skin  (Local  Model) 


1 .  For  each  analysis,  GEOMESH  generates  a  detailed  CAD  model  representing  the  entire  struc¬ 
ture  according  to  the  user  supplied  design  parameters. 

2.  GEOMESH  then  converts  the  CAD  model  into  a  simplified  version  of  the  CAD  model.  The 
simplified  CAD  model  has  the  appropriate  level  of  detail  for  the  current  analysis.  When  an 
intermediate  analysis  is  necessary,  it  can  also  convert  global  analysis  results  to  the  boundary 
conditions  for  the  intermediate  stage.  CAD  models  contain  not  only  geometry  information 
but  also  the  information  required  for  the  finite  element  analysis.  These  parameters  (such  as 
material  properties,  structural  dimensions,  body  forces  and  boundary  conditions)  are  associ¬ 
ated  with  geometric  entities  such  as  vertices,  edges  and  faces. 

3.  Finally,  a  mesh  and  the  associated  boundary  conditions  are  generated  from  the  simplified 
version  of  the  CAD  model. 

To  obtain  boundary  conditions  for  the  local  (or  intermediate  if  necessary)  analysis  from  results 
of  the  global  analysis,  GEOMESH  utilizes  its  post-processing  functionalities  that  extract  analysis 
results  at  any  location  in  the  CAD  model.  Displacements  or  stresses  can  be  obtained  at  any  point 
in  the  model  using  a  finite  element  interpolation  scheme. 

Geometry-Based  Analysis 

Geometry-based  analysis  is  a  relatively  new  concept  in  the  finite  element  analysis  field.  It  allows 
a  user  to  interact  with  a  geometry  model  to  perform  analysis  rather  than  to  manipulate  a  mesh 
directly.  This  is  a  mesh-invisible  approach  wherein  the  mesh  is  hidden  from  the  user.  A  geometry 
model  is  defined  by  a  few  geometrical  parameters.  Associated  analysis  conditions,  such  as  material 
properties  and  boundary  conditions,  are  attached  to  the  geometry  model.  The  analysis  result  can 
be  evaluated  at  an  arbitrary  location  or  region  of  the  geometry  model,  independent  of  the  elements 
and  the  nodes  that  are  hidden  behind  the  geometry  model. 

A  typical  geometry  model  in  an  engineering  field,  especially  the  aircraft  industry,  can  be  de¬ 
fined  using  a  small  number  of  dimensional  and  topological  parameters.  Once  the  geometry  is 
parameterized,  parametric  study  can  be  easily  performed. 

At  first  glance,  geometry-based  analysis  seems  to  have  nothing  to  do  with  hierarchical  anal¬ 
ysis.  However,  subregion  definition  and  simplification  are  highly  problem-dependent  in  a  typical 
hierarchical  analysis.  They  change  frequently  or  even  dynamically  in  run  time  during  a  parametric 
study.  If  the  user  has  to  directly  create  and  manipulate  a  finite  element  model,  the  mesh  genera¬ 
tion  and  boundary  condition  transfer  becomes  extremely  complex  and  time  consuming.  For  this 
reason,  hierarchical  analysis  based  on  traditional  pre-  and  post-processors,  which  require  the  user 
to  manipulate  the  mesh  directly,  is  largely  impractical  for  parametric  study  or  design  optimization. 
Geometry-based  analysis  makes  this  approach  practical. 

Automatic  mesh  generation  and  physical  quantity  evaluation  techniques  are  enabling  technolo¬ 
gies  for  geometry-based  analysis  techniques.  The  technologies  implemented  in  GEOMESH  are 
described  herein. 
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Geometry  Model 

A  geometry  model  is  composed  of  topological  data  and  geometrical  data.  Topological  data  repre¬ 
sents  the  relationship  between  topological  entities  such  as  the  vertex,  edge  and  face.  Geometrical 
data  defines  the  shape  of  topological  entities  using  a  geometrical  entity  such  as  a  surface,  polygon, 
line  and  curve.  A  simplified  form  of  boundary  representation  (B-rep)  is  used  for  the  topological 
data  of  a  geometry  model.  Because  of  its  generality,  the  boundary  representation  has  been  recently 
adopted  in  a  broad  range  of  applications  in  CAD/CAM/CAE  industries. 

A  geometry  entity  defines  the  shape  of  a  topological  entity.  The  straight  line  and  flat  polygon 
are  implemented  as  the  primitive  geometry  entities.  The  shape  of  each  loop  is  defined  by  a  cor¬ 
responding  flat  polygon.  The  mapping  from  a  flat  surface  to  a  cylindrical  surface  is  used  for  the 
modeling  of  a  curved  panel. 

Mesh  Generation 

GEOMESH  can  automatically  generate  meshes  consisting  of  3-D  beam  and  shell  elements  using 
a  mapped  mesh  approach.  The  mapped  mesh  generation  consists  of  two  steps:  i)  generation  of 
super-elements  for  a  geometry  model;  and  ii)  decomposition  of  these  super-elements  into  finite 
elements.  This  approach  generates  a  minimum  number  of  high  quality  quadrilateral  elements  for  a 
relatively  regular  geometry  (which  is  typical  in  a  aircraft  structure). 

In  the  current  implementation,  super-elements  can  be  generated  automatically  only  if  the  ge¬ 
ometry  model  is  composed  of  faces  which  are  flat  and  parallel  to  the  XY,  YZ  or  ZX  planes,  and 
edges  which  are  straight  and  parallel  to  the  X,  Y  or  Z  axis.  This  is  one  of  the  major  limitations  in 
the  current  implementation.  It  will  be  addressed  in  the  future. 

Mapped  Mesh  Approach  In  the  mapped  mesh  approach  adopted  in  GEOMESH,  each  face  of 
a  geometry  model  is  divided  into  super-elements.  The  shape  of  the  super-elements  is  quadrilat¬ 
eral.  A  number  of  divisions  is  specified  in  each  division  direction.  If  the  numbers  of  divisions 
are  inconsistent  among  adjacent  super-elements,  quadrilateral  elements  generated  from  the  map¬ 
ping  become  inconsistent  too.  If  a  crack  is  located  between  adjacent  super-elements,  2  nodes  are 
generated  at  the  same  position  for  the  formation  of  unconnected  elements. 

Automatic  Super-element  Generation  GEOMESH  can  automatically  generate  super-elements 
when  a  geometry  model  is  orthogonal,  (i.e  the  model  is  composed  of  faces,  which  are  flat  and  par¬ 
allel  to  the  XY,  YZ  or  ZX  planes,  and  edges,  which  are  straight  and  parallel  to  the  X,  Y  or  Z  axis.) 
Since  the  global  damage  tolerance  model  satisfies  this  condition,  the  super-element  generation  is 
fully  automated. 

An  algorithm  is  implemented  to  create  a  3-D  orthogonal  non-uniformly  spaced  grid  for  each 
part,  so  that  each  vertex  of  the  part  is  placed  on  one  of  the  grid  points.  Each  face  is  then  classified  as 
to  whether  the  face  is  parallel  to  the  XY,  YZ  or  ZX  plane,  and  a  2-D  cut  plane  which  is  coincident  to 
the  face  is  identified  from  the  3-D  grid.  The  2-D  cut  plane  is  a  collection  of  rectangles.  Rectangular 
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super-elements  which  are  contained  within  the  face  are  collected  from  the  rectangles  of  the  2-D 
cut  planes.  Because  all  the  super  elements  generated  by  the  algorithms  are  rectangles,  it  is  easy  to 
define  the  number  of  divisions.  If  an  average  element  size  is  specified,  the  AGILE  mesh  generator 
calculates  the  number  of  divisions  for  all  the  super  elements.  Then,  it  automatically  generates  shell 
and  beam  elements. 

Evaluation  of  Analysis  Result 

GEOMESH  can  evaluate  analysis  results  at  an  arbitrary  location  independent  of  the  locations  of 
the  nodes  and  elements,  as  if  the  analysis  results  were  mapped  on  the  geometry  model  rather  than 
at  each  node  or  at  each  element  integration  point.  The  displacements  of  the  beams  are  represented 
using  free-form  curves;  and  the  displacements  of  the  shells  are  represented  using  free-form  sur¬ 
faces. 

A  physical  quantity  at  a  given  point  is  evaluated  through  :  i)  searching  for  the  element  that 
contains  the  given  point;  ii)  calculating  the  local  co-ordinate  of  this  point  in  this  element  ;  iii) 
calculating  the  value  of  the  physical  quantity  by  interpolation  either  from  values  at  the  nodes  or  at 
the  element  integration  points  using  standard  finite  element  interpolation  techniques. 

Note  that  this  algorithm  does  not  work  if  an  evaluation  point  is  coincident  with  the  location  of 
a  crack,  because  it  is  ambiguous  which  side  of  the  crack  is  chosen.  Therefore,  the  evaluation  point 
has  to  be  placed  a  small  distance  apart  from  the  crack. 

Displacement  Displacements  are  imported  at  each  node  in  a  finite  element  model.  The  displace¬ 
ments  have  6  components,  which  consists  of  3  translation  components,  Ux,  Uy  and  u^,  and  3  rotation 
components,  Qx,  0);  and  0^.  The  distribution  of  the  displacement  is  assumed  to  be  linear  over  a  4 
node  quadrilateral  shell  element. 

Stress  In-plane  element  stress  is  imported  at  each  integration  point  of  each  shell  element.  The 
local  co-ordinate  system  of  the  element  is  used,  the  in-plane  stress  has  3  components,  (^xy 
Cyy.  The  in-plane  element  stress  is  redistributed  at  each  node  to  form  an  extrapolated  nodal  stress, 
so  that  the  distribution  of  stress  is  smooth.  Currently,  this  is  implemented  in  the  1st  order  4  node 
quadrilateral  shell  element.  This  element  has  one  integration  point  at  the  center  of  the  element. 
Therefore,  the  distribution  of  in-plane  element  stress  is  constant  over  the  element.  However,  the 
distribution  of  extrapolated  nodal  stress  is  linear  over  the  element. 

Reaction  Force  GEOMESH  uses  beam  elements  to  model  stiffeners,  rivets  or  rigid  connections. 
A  rivet  force  may  be  needed  at  the  point  where  a  rivet  is  directly  connected  to  a  skin  to  form  free 
body  diagrams. 

A  reaction  force  is  calculated  from  the  corresponding  beam  elements.  For  each  beam  element, 
a  shear  force  is  calculated  from  the  displacements  of  the  2  nodes  of  the  beam  element. 
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233  Analysis  Model  and  Feature  Modeling 

An  analysis  model  is  a  geometry  model  with  associated  attributes  relevant  to  a  finite  element  anal¬ 
ysis.  Known  as  a  Product  Modefi  in  the  CAD  industry,  it  is  an  intelligent  CAD  model  with  infor¬ 
mation  and  knowledge  about  the  finite  element  analysis.  An  analysis  model  is  a  product  model  for 
the  finite  element  analysis.  It  is  the  primary  modular  interface  for  the  finite  element  analysis. 

In  general,  the  user  can  create  an  analysis  model;  interact  with  it;  attach  analysis  conditions  to 
it;  and  visualize  analysis  results  on  it.  Because  it  is  completely  hidden  behind  the  geometry  model, 
the  mesh  is  invisible  from  the  user. 

A  feature  model  is  a  high  level  geometry  model  with  built-in  information  of  model  simplifica¬ 
tion.  It  consists  of  a  history  list  of  geometry  definition  commands.  Each  command  is  an  operation 
such  as  extrusion,  fillet,  cut-out,  or  so.  A  geometry  model  of  a  specified  level  of  detail  can  be 
constructed  from  its  command  history  list. 

A  feature-based  model  knows  how  to  analyze  itself  and  how  to  perform  model  simplification 
if  necessary.  It  can  generate  an  analysis  model.  It  also  knows  how  to  transfer  boundary  conditions 
between  analysis  models  of  different  levels  of  details. 

The  part  and  connection  are  the  building  blocks  of  a  feature  model.  A  feature  model  consists 
of  parts  and  connections.  A  connection  joins  two  parts  to  each  other.  Each  part  and  connection 
contains  a  command  history  list  to  generate  a  geometry  model  of  the  most  detailed  level.  By 
specifying  the  level  of  detail  on  each  part  connection,  the  feature  model  can  produce  a  geometry 
model  that  is  simplified  to  the  given  level  of  detail. 

Part  Two  types  of  parts  for  feature  models  are  implemented.  They  are:  beam-like  narrow  parts 
and  shell-like  thin  parts. 

1.  Beam-like  narrow  parts  can  produce  a  wire  frame  model  of  a  beam  structure  or  a  surface 
model  of  a  shell  structure.  The  knowledge  to  transfer  boundary  conditions  between  two  wire 
frame  models,  a  wire  frame  model  and  a  surface  model,  or  two  surface  models  has  been 
implemented. 

2.  Shell-like  thin  parts  can  produce  a  surface  model  of  shell  structure.  The  knowledges  to 
transfer  boundary  conditions  between  two  surface  models  are  implemented. 

Connection  Three  types  of  connections  are  implemented.  They  are:  riveted  line,  riveted  area, 
and  rigid  connection.  In  a  riveted  line,  parts  are  connected  by  equally  spaced  rivets.  In  a  riveted 
area,  parts  are  connected  by  a  rectangular  grid  of  rivets.  The  rectangle  is  perpendicular  to  the 
direction  of  the  rivets. 

^  A  product  model  is  a  CAD  geometry  model  which  contains  application-specific  information  and  intelligence  to 
perform  application-specific  tasks. 
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2.3.4  Hierarchical  Analysis  Based  on  Feature  Modeling 

In  an  automated  hierarchical  analysis,  only  a  single  detailed  feature  model  is  required.  The  hi¬ 
erarchical  analyzer  automatically  constructs  the  analysis  models  for  each  level  of  analysis  with 
corresponding  geometrical  details  and  coverage.  Thus,  four  major  tasks  are  performed  by  the 
automated  hierarchical  analyzer.  They  are  i)  the  simplification  of  geometry;  ii)  the  extraction  of 
subregions;  iii)  the  evaluation  of  the  analysis  results  in  the  current  analysis;  and  iv)  the  construction 
of  boundary  conditions  for  the  analysis  in  the  next  stage.  These  tasks  are  described  in  this  section. 

Simplification  of  Geometry 

Simplification  of  geometry  is  necessary  to  reduce  the  number  of  D.O.F.  required  to  construct  a 
finite  element  model  for  a  given  feature  model.  The  D.O.F.' s  saved  by  simplifying  unnecessary 
details  can  be  allocated  to  the  regions  where  mesh  refinement  can  significantly  affect  the  accuracy 
of  the  finite  element  analysis,  such  as  the  vicinity  of  a  crack  tip. 

In  general,  simplification  of  geometry  is  the  process  of  elimination  of  small  scale  dimensional 
parameters  in  the  feature  model.  Small  features  can  be  eliminated  without  significant  degradation 
in  accuracy  (except  in  the  vicinity  of  the  eliminated  feature).  For  example,  in  case  of  an  aircraft 
fuselage,  dimensional  parameters  such  as  the  radius  of  rivet  holes  in  the  skin,  the  size  of  MSD 
cracks,  a  joggle  depth  of  stringers  at  the  junction  with  tear  straps,  and  the  radius  of  fillets  of  a 
finger  doubler  can  be  ignored  in  a  typical  intermediate  level  model. 

Narrow  structures  can  be  simplified  as  beams.  If  all  the  stiffeners  are  modeled  as  beams  in 
a  global  model,  the  dimensional  parameters  in  the  same  order  or  less  than  that  of  the  width  (or 
height)  of  the  stiffeners  are  ignored  from  the  feature  model.  The  information  is  used  to  generate 
the  geometrical  properties  of  the  beams. 

The  rivet  distribution  pattern  heavily  affects  the  number  of  D.O.F.  required  to  model  a  finite 
element  model.  To  reduce  the  number  of  D.O.F.,  rivet  spacing  can  be  artificially  enlarged.  At  the 
same  time,  the  geometrical  property  parameters  of  the  rivets  have  to  be  modified  accordingly. 

Model  simplification  methods  that  have  been  implemented  are  i)  modeling  stiffeners  as  beams; 
ii)  elimination  of  small  features;  iii)  modeling  rivets  as  beams;  iv)  simplification  of  rivet  distribu¬ 
tions;  v)  using  rigid  connections. 

Stiffener  Modeled  as  a  Beam  A  narrow  structure  is  defined  by  an  extrusion  operation  in  CAD 
geometry  modeling.  First,  a  geometry  model  of  a  cross-section  is  defined.  Then,  it  is  extruded 
along  a  path  which  defines  the  shape  of  the  narrow  structure.  A  solid  model  is  generated  from  the 
extrusion  of  a  surface;  a  surface  model  is  generated  from  the  extrusion  of  a  wire  frame;  and  a  wire 
frame  is  generated  from  the  extrusion  of  a  point. 

A  stiffener  is  the  extrusion  of  a  section  along  a  given  path.  It  can  be  modeled  as  a  3D  solid 
with  a  solid  model,  or  as  a  shell  structure  with  a  surface  model,  or  as  a  frame  with  a  wire  frame 
model.  In  order  to  simplify  the  stiffener  as  a  shell  structure,  the  thickness  parameter  is  degenerated 
into  a  geometrical  property  of  a  surface  model.  In  order  to  simplify  the  stiffener  as  a  frame,  the 
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cross-sectional  parameter  is  degenerated  into  a  geometrical  property  of  a  wire  model.  A  feature 
model  of  the  stiffener  is  implemented  so  that  the  automatic  hierarchical  analyzer  can  perform 
these  model  simplifications.  Automatic  mesh  generators  have  been  implemented  to  generate  shell 
models  and/or  beam  models.  Since  no  solid  mesh  generator  is  implemented,  no  3D  solid  model 
can  be  generated. 

A  wire  frame  model  is  used  for  a  stiffener  in  a  aircraft  structure  in  a  global  analysis.  It  is 
discretized  using  beam  elements  in  the  corresponding  finite  element  model. 

Elimination  of  a  Small  Feature  The  complexity  of  a  geometry  model  affects  significantly  the 
number  of  D.O.F.  of  the  generated  finite  element  model.  Some  of  the  small  features  can  be  ignored 
in  a  early  stage  of  a  hierarchical  analysis,  since  the  detailed  stress  solution  nearby  these  features  is 
not  important.  In  the  early  stage  of  a  hierarchical  analysis,  it  is  important  to  determine  the  accurate 
boundary  conditions  for  the  analysis  in  the  next  stage.  These  small  features  have  an  insignificant 
effect  on  the  solution  for  the  boundary  condition;  while  accurate  modeling  of  these  features  can 
significantly  increase  the  number  of  D.O.F. 

Filleting,  tapering,  making  a  small  hole,  or  offsetting  are  examples  of  such  small  features.  Even 
small  curved  edges  can  be  considered  as  features,  because  the  extra  nodes  needed  to  represent 
the  curve  can  lead  to  a  significant  increase  in  the  number  of  D.O.F.  in  the  finite  element  model 
generated  by  the  automated  mesh  generator. 

Typically,  the  construction  of  a  geometry  model  starts  from  a  simple  model.  Smaller  and 
smaller  features  are  added  incrementally  to  it.  With  a  record  of  this  procedure,  a  simplified  model 
can  be  reproduced  by  a  modification  to  the  procedure. 


Rivet  Modeling  as  a  Beam  Beam  elements  are  used  to  model  the  rivets.  Swift's  experimental 
equation  is  used  to  determine  the  stiffness  of  the  rivets.  They  are 
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where  Eaxiai  is  the  axial  direction  stiffness  of  the  rivet  {psi  ■  in),  Etrans  is  the  transverse  direction 
stiffness  of  the  rivet  (psi  •  in),  E  is  the  Young's  modulus  of  the  rivet  (psi),  D  is  the  diameter  of  the 
rivet  (in),  tsi  is  the  thickness  of  sheet  1  (in),  and  ts2  is  the  thickness  of  sheet  2  (in).  A  and  C  are  the 
coefficients  in  Swift's  equation.  A  =  5.0,  C  =  0.8  for  aluminum  rivets;  and  A  =  1.666,  C  =  0.86  for 
steel  rivets. 

The  corresponding  geometrical  properties  for  the  beam  elements  are  defined  as  follows. 
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where  A  is  the  cross  sectional  area  (in^),  I  is  the  moment  of  inertia  and  /  is  the  length  of  beam 

element  of  the  rivet  (in). 

Simplification  of  the  Rivet  Distribution  A  feature  model  of  a  rivet  connection  consists  of  a 
number  of  rows  of  rivets.  These  rows  of  rivets  are  distributed  in  a  regular  fashion.  Multiple  rivets 
can  be  combined  and  modeled  by  a  single  beam  element,  whose  stiffness  is  the  sum  of  these 
rivets.  Multiple  rows  of  rivets  can  be  combined  and  modeled  by  a  single  row  of  beams.  This 
simplification  technique  artificially  enlarges  the  rivet  spacing.  Therefore,  the  number  of  D.O.F.  in 
the  finite  element  model  can  be  reduced  significantly.  This  is  a  very  useful  simplification  technique 
for  a  earlier  stage  analysis  of  a  build-up  structure  in  a  hierarchical  analysis.  In  an  extreme  case,  all 
the  rivets  in  the  feature  model  can  be  degenerated  to  a  single  beam  element  located  at  the  center  of 
the  rectangle. 

Rigid  Connection  as  a  Collection  of  Rivets  The  connection  between  a  frame  and  a  stringer, 
without  a  shear  clip,  is  effectively  a  rigid  connection.  The  geometrical  details  of  such  a  connection 
are  very  complicated.  Such  a  connection  is  modeled  as  a  rigid  connection,  using  very  strong  beams. 
The  stiffness  of  the  beam  element  is  assigned  by  the  program  so  that  the  beam  element  is  so  stiff 
that  the  relative  translation  and  rotation  between  the  two  nodes  are  practically  zero. 

Subregion  Extraction 

A  subregion  is  automatically  extracted  for  the  subsequent  models  involved  in  a  hierarchical  anal¬ 
ysis.  A  subregion  can  be  of  any  shape.  Typically  it  is  rectangular.  Usually,  the  subregion  is 
contained  completely  inside  the  given  model,  from  which  the  extraction  is  performed.  However, 
the  prescribed  subregion  can  also  intersect  with  the  boundary  of  the  given  model.  In  this  case, 
the  boundary  of  the  given  model  in  the  prescribed  subregion  becomes  a  part  of  the  boundary  in 
the  extracted  model.  A  CAD  operation,  intersection  operation,  is  used  to  extract  the  geometrical 
model  of  the  subregion. 

Construction  of  Boundary  Conditions 

The  boundary  conditions  for  the  subsequent  analysis  are  constructed  using  the  analysis  results 
evaluated  at  the  locations  of  the  boundaries  in  the  current  model.  Once  the  subregion  extraction  is 
performed,  the  location  of  the  boundaries  in  the  current  model  is  determined.  Analysis  results  can 
be  evaluated  at  these  locations  for  the  calculation  of  the  boundary  conditions  for  the  analysis  in  the 
next  stage. 

The  boundary  conditions  at  the  end  of  stiffeners  are  prescribed  as  displacement  boundary  con¬ 
ditions  in  the  subsequent  models  generated  by  the  automatic  hierarchical  analyzer.  When  a  wire 
frame  model  of  a  stiffener  is  modeled  in  detail  using  a  surface  model  in  the  next  stage,  the  rotation 
and  the  displacement  at  this  point  is  used  to  calculate  the  displacement  constraints  of  the  shell 
model  at  the  end  of  the  stiffener. 
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When  the  hierarchical  analysis  transitions  from  a  3D  shell  analysis  to  a  2D  plane  stress/strain 
analysis,  it  is  necessary  to  perform  a  projection  of  the  analysis  result.  For  example,  when  a  local 
region  of  the  fuselage  skin  is  to  be  analyzed  by  a  2D  plane  stress/strain  analyzer  (such  as  a  code 
using  a  2D  finite  element  alternating  method),  conversion  of  analysis  results  for  a  3-D  shell  to  the 
corresponding  boundary  condition  for  a  2D  analysis  is  required.  In-plane  displacements  in  the  mid¬ 
plane  are  used  for  displacement  boundary  conditions;  while  membrane  stresses  in  the  mid-plane 
are  used  for  the  construction  of  traction  boundary  conditions. 

2.4  Buckling  Analysis  of  a  Composite  with/without  deiamination 

The  primary  damage  in  a  structure  loaded  in  compression  which  must  be  considered  during  design 
is  a  hole  (Fig.  2.12).  For  laminated  composites,  delaminations  will  also  be  of  relevance.  It  is 
to  be  expected  [with  some  preliminary  results  that  confirm  these  expectations  being  reported  in 
Vellaichamy,  et.  al  (1990)]  that  the  buckling  load  of  a  panel,  with  a  hole,  in  compression  will 
depend  on  the  shape  and  the  size  of  the  hole.  The  buckling  load  of  a  panel  with  an  elliptical  hole 
whose  major  axis  is  aligned  with  the  direction  of  compression  can  be  expected  to  be  lower  than 
that  with  a  circular  hole  of  the  same  area.  Further,  it  is  expected  that  for  an  elliptical  hole  of  a  given 
area,  the  buckling  load  will  decrease  with  an  increase  in  the  aspect  ratio  of  the  ellipse,  as  long  as  the 
major  axis  of  the  ellipse  is  aligned  with  the  direction  of  compression.  In  the  limit  as  the  elliptical 
hole  shrinks  to  a  crack  whose  axis  is  parallel  to  the  direction  of  compression,  with  the  slightest 
imperfection,  the  crack  may  propagate  in  mode  III  in  post-buckling  deformation.  Such  mode 
III  behavior  would  severely  affect  the  structural  integrity.  BUCKDEL,  as  developed  primarily  in 
Phase  I  of  this  SBIR  project,  can  be  used  to  compute  the  buckling  and  post-buckling  response  of 
a  stiffened  structure  containing  damage  in  the  form  of  holes  and  delaminations  of  arbitrary  shape. 
In  addition,  it  computes  the  pointwise  energy  release  rate  around  the  delamination  front.  In  future 
work,  the  ability  to  treat  the  mode  III  crack  problem  will  be  added  to  BUCKDEL. 

A  brief  overview  of  BUCKDEL  is  given  here.  For  more  details,  see  the  Users  and  Theory 
Manuals  of  BUCKDEL,  which  are  parts  of  the  final  report  for  Phase  I  of  this  project.  BUCKDEL 
was  developed  as  a  stand-along  program  in  Phase  I  of  this  project.  It  was  connected  with  ASTROS 
as  a  linear  buckling  analysis  module  in  Phase  II. 

BUCKDEL  uses  a  multi-domain  method  to  model  delaminations  of  arbitrary  shape.  In  this 
method,  the  delaminated  shell  is  assumed  to  be  assembled  with  three  regions-(l)  Undelaminate: 
undelaminated  zone;  (2)  Delaminate:  thinner  side  of  the  delaminated  zone  and  (3)  Base:  thicker 
side  of  the  delaminated  zone.  Transverse  shear  deformation  plays  an  important  part  in  the  case  of 
composite  structures,  hence,  it  is  introduced  explicitly  and  the  assumptions  of  the  Reissner-Mindlin 
theories  of  plate  bending  are  used  for  modeling  each  region  of  the  multi-domain  model.  Thus,  for 
each  region,  the  3-dimensional  displacement  field  (U  =  {Ui  U3})  is  expressed  in  terms  of  the 
corresponding  mid-surface  displacement  (u  =  {wi  M2  ^3})  and  rotation  (6  =  {0i  62  0})  fields  as, 

uW(xa,X3)  =  u(')(xa)  -4'^0^'^(xa)  (2.5) 

where  (a  =  1 , 2)  are  the  inplane  curvilinear  shell  coordinates  and  X3  ^  is  the  thickness  coordinate 
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Figure  2.12:  Damage  in  the  Upper  (Compression)  Skin  of  a  Wing 


for  the  {i  =  1,2,3)  shell.  The  structural  continuity  at  the  delamination  front  F  is  maintained  by 
assuming  the  deformation  to  be  unique  at  the  junction  of  the  three  shells  i.e. 
on  r  in  accordance  with  the  Reissner-Mindlin  law  of  flexure  (Eq.  (2.5)).  In  other  words,  at  the 
delamination  edge,  the  mid-surface  degrees  of  freedom  of  the  delaminate  and  the  base  shells  are 
assumed  to  be  related  to  those  of  the  undelaminated  shell  by, 

eS,”  =  =  e«  (2.6) 

\  «!;)  =  /„r 

where  is  the  distance  of  the  midsurface  of  the  i‘^  shell  from  the  laminate  midsurface. 

Each  lamina  is  assumed  to  be  orthotropic,  and  the  inplane  stresses  =  {gh  G22  ori2}^*^  and 
the  transverse  shear  stresses  =  {ti3  X23}^'^  are  related  to  the  linear  components  of  the  membrane 

strain  =  {sn  £22  £12}^'^  nonlinear  components  of  the  membrane  strain  =  {vu  V22 
flexural  strain  due  to  mid-surface  rotation  k(')  =  {kh  K22  Ki2}('h  flexural  strain  due  to  transverse 
shear  strain  =  {xn  X22  Xi2}^'^  and  transverse  shear  strains  =  {Yi3  723}^'^  as 

£^11  ^12  ^16  0  0 

£"12  E22  E26  0  0 

£^16  E26  £^66  0  0 

0  0  0  £44  £45 

0  0  0  £45  £55 

where  the  material  constitutive  terms  E^‘)  are  functions  of  the  thickness  coordinate  of  each  shell 

X3  Generally,  for  a  laminate  with  orthotropic  layers,  E^j)  are  assumed  to  be  piecewise  constants 
over  the  laminate  thickness. 

Integrating  along  the  thickness,  the  constitutive  equations  can  be  written  in  terms  of  the  inplane 
stress  resultants  N  =  {A^ii  N22  Nn}^  bending  moments  M  =  {Mu  M22  M12}.  and  transverse  shear 
stress  resultant  Q  =  {<2i3  <223 }>  fo**  sach  region  of  the  multi-plate  model  as, 

M  i  =  B  D  0  i  (K+x)  >  (2-8) 

qJ  Lo  0  gJ  i  Y  > 

where 

(Aki-,  Bkf,  Dkif^  =  J^Elf{x3){l;x3;4)dx3 

=  f  SmSnE\f{X3)dx3 

In  addition  to  a  beam  element,  BUCKDEL  implements  a  three  noded  triangular  curved  shell 
element.  The  shell  element  is  described  in  the  curvilinear  coordinate  system  x  —  y  and  the  area 
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coordinates  are  used  for  the  field-description.  Accordingly  we  have, 

3 

=  1}/  (2.9) 

/=l 

Inverting  the  above  relationship,  we  get, 

Li  =  ~  {anx  +  any  +  )  (2. 10) 

where 

an  =  yj  -  yk\  aa  =Xk-  xj;  uq  =  xjy^  -  x^yj 
^  ^  (^2}'3  -  ^3>'2  +  x^yi  -Xiy2+  Xiy2  -  X2yi ) 

and  7  =  2, 3, 1 ;  ^  =  3, 1 , 2  as  /  =  1 , 2, 3. 

The  inplane  displacements  and  the  transverse  shear  strains  need  to  satisfy  C®-continuity  while 
the  transverse  deflection  needs  to  satisfy  -continuity  in  the  present  formulation.  The  independent 
field  variables  u,  v,  w,  yxz  and  are  expressed  in  terms  of  the  nodal  degrees  of  freedom  Ui,  Vj,  w,, 
Pxi  =  (-w,y  )i,  py,.  =  {w,x  )i,  Jxzi  and  jy^.  as. 


3 

{uvyxzYyz}  =  ^Li{uvyxzyyz}i 
3 

w  =  £(iVi,w/  +  A2,K+iV3/p).,)  (2.11) 

where 


are  the  cubic  polynomials  for  the  transverse  deflection. 

In  the  above  element  formulations,  the  inter-element  C®-continuity  is  exactly  satisfied  for  all 
the  field  variables.  However,  the  inter-element  C’  -continuity  required  for  the  transverse  deflec¬ 
tion,  in  case  of  a  shell  element,  is  satisfied  a  posteriori  in  a  weak  form  using  the  Hu-Washizu 
variational  principle. 

BUCKDEL  uses  an  automated  method  to  follow  the  post-buckling  paths  of  the  damaged  struc¬ 
ture.  Teh  automated  post-buckling  solution  involves:  detection  of  a  possible  instability  in  the 
solution  and  elimination  of  possible  path-retracing;  classification  of  the  detected  instabilities;  and 
computation  of  the  post-through  buckling  solution(s).  Solution  instabilities  are  detected  by  moni¬ 
toring  the  rank  of  the  tangent  stiffness  matrix.  Whenever  the  determinant  of  the  tangent  stiffness 
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matrix  changes  its  sign,  the  solution  senses  possible  instabilities  in  that  range  of  load  and  changes 
the  sign  of  the  next  load  increment  to  avoid  path-retracing.  Through  a  cycle  of  iterations,  the 
location  of  instabilities  are  identified  as  the  load  levels  for  which  the  tangent  stiffness  becomes  sin¬ 
gular.  The  tangent  stiffness  is  often  scaled  to  minimize  numerical  errors.  The  identified  instability 
points  are  then  classified  as  limit  points  or  bifurcation  points  using  some  simple  and  cost-effective 
rules  [Huang  and  Atluri  (1995)].  If  the  instability  point  is  a  limit  point,  the  arc-length  continua¬ 
tion  is  enough  to  obtain  the  post-buckling  solution  path.  However,  if  the  instability  point  happens 
to  be  a  bifurcation  point,  the  strategies  described  in  detail  in  Huang  and  Atluri  (1995)  are  used 
to  trace  the  appropriate  post-buckling  solution  branch.  The  nonlinear  fundamental  state  between 
the  two  solution  points  n  -  1  and  n  in  the  neighborhood  of  the  bifurcation  point  is  linearized  to 
obtain  the  asymptotic  solution  for  obtaining  an  approximate  critical  buckling  load  factor.  A  linear 
combination  of  the  normalized  eigen-vector  associated  with  the  critical  buckling  load  factor  and 
its  orthogonal  counterpart  is  used  to  determine  the  initial  post-buckling  paths. 

2.5  Finite  Element  Alternating  Method 

The  Schwartz-Neumann  alternating  method  is  based  on  the  superposition  principle.  The  solution 
on  a  given  domain  is  the  sum  of  the  solutions  on  two  other  overlapping  domains,  with  part  of 
the  boundary  conditions  as  unknowns.  The  alternating  method  can  be  viewed  as  the  fixed  point 
iteration  scheme  used  to  solve  these  unknown  boundary  conditions.  Based  on  this  point  of  view,  we 
can  perform  a  convergence  study.  The  alternating  method  converges  unconditionally  when  there 
are  only  traction  boundary  conditions  specified  on  the  body.  The  convergence  criterion  for  mixed 
boundary  value  problems,  where  there  are  applied  displacement  boundary  conditions  as  well  as 
traction  boundary  conditions,  is  discussed  in  the  following.  Compare  the  work  done  by  the  applied 
forces  in  the  following  two  cases.  In  the  first  case,  arbitrary  displacement  conditions  exist  on  the 
surfaces  of  the  cracks  in  the  cracked  finite  body,  while  all  the  boundary  conditions  elsewhere  are 
replaced  by  homogeneous  boundary  conditions,  i.e.  remove  all  tractions  and  reduce  all  the  applied 
displacements  to  zero  magnitude.  In  the  second  case,  the  same  displacement  conditions  exist  on  the 
surfaces  of  the  cracks  in  the  infinite  domain.  If  the  work  done  in  the  cracked  finite  body  is  always 
smaller  than  twice  the  work  done  in  the  infinite  domain,  then  the  alternating  method  converges. 
Otherwise,  it  does  not.  For  most  practical  problems,  this  ratio  is  close  to  one.  Thus,  the  alternating 
method  converges  rapidly,  as  discussed  in  detail  in  the  following  section. 

2.5.1  Superposition  Principle  and  the  Alternating  Method 

Consider  n  cracks  in  a  body  of  finite  size.  The  crack  surfaces  which  are  traction  free,  are  de¬ 
noted  collectively  as  Tc.  Let  the  boundary  of  the  finite  domain(not  including  the  crack  surface) 
be  r,  of  which  the  boundary  with  prescribed  tractions  t"  is  Ff,  and  the  boundary  with  prescribed 
displacements  u°  is  r„.  It  is  clear  that  F  =  r„  U  Ff . 

The  alternating  method  uses  the  following  two  simpler  problems  to  solve  the  original  one.  The 
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Figure  2.13:  Superposition  Principle  for  Finite  Element  Alternating  Method 

first  one,  denoted  as  Pana  [shown^  in  Fig.  2.13(c)],  is  that  of  the  same  n  cracks  in  the  infinite  do¬ 
main  subjected  to  the  unknown  crack  surface  loading  T.  The  second  one,  denoted  as  Pfem  [shown 
in  Fig.  2.13(b)],  has  the  same  finite  geometry  as  in  the  original  problem  except  that  the  cracks  are 
ignored.  The  boundary  r„  of  Pfem  h^s  the  prescribed  displacement  u,  while  the  boundary  F;  has 
the  prescribed  traction  t.  The  prescribed  displacements  and  tractions  are  different  from  those  in 
the  original  problem  in  general.  Because  of  the  absence  of  the  cracks,  the  problem  Pfem  can  be 
solved  much  easier  by  the  finite  element  method  (or  the  boundary  element  method). 

To  solve  the  original  problem,  Pqrg  (shown  in  Fig.  2.13(a)),  the  crack  surface  loading  T,  the 
prescribed  displacement  u  and  the  traction  t  must  be  found  such  that  the  superposition  of  the  two 
alternative  problems  Pana  and  Pfem  yields  the  original  one,  Pqrg-  The  detailed  procedures  to  find 
these  boundary  conditions  are  described  in  the  following. 

In  the  uncracked  body  problem  Pfem^  the  tractions  T  at  the  location  of  the  cracks  in  the  cracked 
body  Pqrg  can  be  solved,  for  any  given  boundary  loads  u  and  t,  using  the  finite  element  method. 
Due  to  the  linearity  of  the  problem,  the  solution  can  be  denoted  as 

T  =  K“u  +  K‘t  (2.13) 

where  K“  and  K‘  are  linear  operators. 

^  Fig.  2.13  only  illustrates  one  crack.  Many  cracks  may  be  present. 
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Similarly,  the  tractions  on  boundary  F,  and  the  displacements  on  boundary  r„  can  be 
found  in  the  infinite  domain  Pana  for  the  given  crack  surface  load  T,  which  is  the  same  as  the 
crack  surface  traction  obtained  in  the  Pfem-  The  solution  can  be  denoted  as 


if  =  K'^T 


(2.14) 


f  =  K'T 


(2.15) 


where  ^  and  are  also  linear  operators.  Subtract  the  solution  for  Pana  from  the  one  for  Pfem- 
The  resulting  solution  has  zero  tractions  at  the  location  of  the  crack  surfaces.  To  ensure  that  the 
resulting  solution  has  the  same  boundary  conditions  on  F,  the  Eq.  (2.16)  and  Eq.  (2.17)  must  be 
satisfied. 


u  =  f  +  iP  (2.16) 

t=:f  +  f  (2.17) 


The  unknown  tractions  T,  t  and  unknown  displacement  u  can  be  solved  using  these  equations 
[Eq.(2.13)  through  Eq.  (2.17)].  Eliminate  m,  if  and  t,  f  by  substituting  Eq.(2.16),  Eq.(2.17), 
Eq.  (2.14)  and  Eq.  (2.15)  into  Eq.  (2.13)  to  obtain  the  following  equation  for  the  traction  T. 


I-  (k"K"  +  K^K‘'J 


T=  + 


(2.18) 


Eliminate  f  and  T  to 
displacement  u. 


obtain  the  following  equation  for  the  unknown  traction  t  and  unknown 
(/-A){  (2.19) 


where 


K^K'*  K‘*K‘ 

¥k“  Wk‘ 


and  I  is  the  identity  operator. 

Similarly,  we  can  obtain  the  following  linear  system  for  traction  f  and  displacement 


{I-A)X  =  Y 


(2.20) 


where 


Y  = 


¥k^‘  ¥k‘  I  J 


It  is  possible  to  solve  these  equations  directly  to  obtain  the  tractions  T,  t  and  displacement  u. 
But  this  involves  the  evaluation  of  K'*  and  K\  which  requires  solving  the  traction  T  at  the  location 
of  the  uncracked  body  subjected  to  all  different  loading  patterns  u  and  t.  We  have  to  solve  the 
uncracked  body  problem  a  larger  number  of  times,  of  the  same  order  as  that  of  the  total  number  of 
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degrees  of  freedom  of  the  boundary  nodes,  using  the  finite  element  method.  Thus,  it  can  be  very 
expensive  to  find  X  by  solving  directly  the  linear  system  {I  —  A)X  =  T.  A  fixed  point  iteration 
scheme  can  be  used  to  solve  this  linear  system.  The  iterative  scheme  can  be  devised  as: 

/  =  0,l,2,...,oo  (2.21) 

where  .  If  this  procedure  converges,  the  solution  is 


X  = 


1=1 


Since 

the  iterative  scheme  Eq.  (2.21)  is  equivalent  to  the  following  alternating  scheme 

7’W  =  /i:“aW+is:M')  (2.22) 


{ 


(2.23) 


for  /  =  0, 1,2, . .  .,00.  In  this  case,  the  uncracked  body  problem  is  solved  only  a  few  times,  because 
this  fixed  point  iteration  scheme  converges  quickly  for  practical  problems.  Therefore,  the  alternat¬ 
ing  method  is  much  more  efficient  than  solving  the  linear  system  directly.  But  it  should  be  noticed 
that  it  may  not  be  necessary  to  use  the  alternating  method  in  some  cases.  It  can  be  more  efficient 
and  accurate  to  solve  directly  when  multiple  crack  solutions  are  constructed  from  that  for  a  single 
crack.  This  will  be  discussed  in  detail  in  a  later  section. 


2.5.2  Convergence  of  the  Alternating  Method 

First  it  is  shown  that  7  —  A  is  not  singular  and  the  linear  system  Eq.  (2.20)  has  a  unique  solution. 
Suppose  7  -  A  is  singular.  Then,  there  must  exist  a  non-zero  X  such  that  (7  -  A)X  =  0,  which  means 
that  there  exist  non-zero  and  t^,  and  therefore  a  non-zero  T,  such  that 

T  =  K‘*u‘‘-\-K‘f 

u‘^  =  K^T 
f  =  K‘T 

In  this  case  the  analytical  solution  and  the  finite  element  solution  have  the  same  displacement  m" 
on  boundary  r«  and  the  same  traction  r“  on  boundary  T;.  Subtracting  the  analytical  solution  from 
the  FEM  solution,  we  obtain  the  solution  for  the  following  problem.  The  entire  boundary  F  is  free 
of  external  loadings  as  well  as  the  crack  surfaces.  But,  the  FEM  solution  gives  zero  displacements 
for  the  crack  surfaces,  while  the  analytical  solution  gives  non-zero  displacements  for  the  crack 
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Figure  2.14:  Subtract  XPfem  from  to  Obtain  the  Solution  for  Pres  which  has  Homogeneous 
Boundary  Condition  on  F 


surfaces  because  of  the  non-zero  T.  Thus,  the  resulting  solution  has  non-zero  displacements  at  the 
crack  surfaces.  This  is  a  contradiction  because  the  cracks  cannot  be  opened  without  any  external 
load.  Consequently,  /  —  A  is  not  singular. 

The  fixed  point  iteration  scheme  Eq.  (2.21)  converges  if  all  the  eigenvalues  of  A  are  in  the  open 
interval  (-1,1).  The  scheme  of  Eq.  (2.21)  converges  since  the  eigenvalues  of  A  are  in  (-1, 1) 
for  most  problems  of  practical  interest.  The  eigenvalues  of  A  are  smaller  than  1 .  Let  Xx  be  an 
eigenvector  of  A  corresponding  to  the  eigenvalue  X. 

T  =  K"{u^)  +  K\t'^) 

Xu^  =  K^T 
Xt^  =  ¥T 

The  solution  Pres,  shown  in  Fig.  2.14(a),  is  obtained  by  subtracting  X  times  the  FEM  solu- 
tion(Fig.  2.14(c))  from  the  analytical  solution(Fig.  2.14(b)).  Here,  u  =  0  and  t  =  0  on  F  and 
the  crack  surface  loading  is  (1  -  X,)r,  while  the  displacements  at  the  crack  surfaces  are  the  same 
as  those  in  the  analytical  solution.  If  the  work  done  in  opening  the  cracks  in  the  infinite  domain 
is  W,  the  work  done  in  opening  the  cracks  in  the  finite  domain(with  the  boundary  condition  u  —  0 
and  t  =  0)  is  ( 1  -  X)1F,  which  is  equal  to  the  strain  energy  stored  in  the  body.  It  must  be  positive. 
Thus,  X.  <  1. 
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It  can  be  shown  that  A,  >  0  in  the  absence  of  the  prescribed  displacement  boundary  conditions. 
In  this  case,  the  resulting  solution  from  the  subtraction  has  zero  tractions  at  the  boundary  T.  Apply 
additional  load  h  to  the  boundary  T  with  the  crack  surfaces  fixed.  The  stress  state  in  the  body 
will  be  the  same  as  that  in  the  analytical  solution  described  above,  after  this  additional  loading  is 
applied.  This  procedure  of  adding  load  on  the  boundary  T  is  exactly  the  same  as  that  in  the  FEM 
solution  for  the  uncracked  body,  except  that  the  load  level  is  X  times  that  in  the  FEM  solution, 
because  the  crack  surfaces  are  fixed.  Therefore,  the  work  done  by  the  additional  load  is  positive. 
Consequently,  ( 1  —  A.)  IT  <  IT  and  A,  >  0.  So,  the  alternating  method  converges  for  cracks  in  finite 
domains  with  arbitrary  shapes  and  arbitrary  traction  boundary  conditions. 

In  general,  the  eigenvalue  A,  can  be  smaller  than  zero  for  mixed  boundary  problems.  It  is 
greater  than  -1  only  if  ( 1  —  A.)  IT  <  2IT.  Thus,  the  convergence  criterion  for  the  alternating  method 
for  the  general  case  with  mixed  boundary  conditions  can  be  stated  as  follows.  The  alternating 
method  [Eq.  (2.21)]  converges  if  the  crack  surface  loads  do  less  work  in  the  finite  domain,  with 
the  homogeneous  boundary  condition  m  =  0  and  t  =  0  on  F,  than  twice  as  much  as  they  do  in  the 
infinite  domain  for  any  arbitrary  distribution  of  crack  surface  displacements  (see  Fig.  2.15). 

Quick  convergence  can  be  expected  for  most  of  the  practical  applications.  For  any  crack  sur¬ 
face  displacements,  the  displacements  and  stresses  at  a  point  decay  rapidly  as  the  point  moves 
away  from  the  cracks.  Thus,  the  work  done  in  the  finite  domain  with  the  homogeneous  boundary 
condition  is  very  close  to  the  work  done  in  the  infinite  domain.  This  implies  that  the  eigenvalues 
of  A  are  very  small  and  the  fixed  point  iteration  converges  rapidly.  Indeed,  all  the  mixed  boundary 
value  problems  we  have  solved  (for  both  2D  and  3D  problems)  to  date  using  the  finite  element 
alternating  method  have  converged. 

2.5.3  Summary  of  FEAM  Procedure 

The  alternating  procedure  defined  in  Eq.  (2.21)  can  be  translated  into  the  following  simple  proce¬ 
dure.  Refer  to  Fig.  2.13. 

1.  Solve  PpEM  with  the  given  load  on  the  boundary  F.  Solve  for  the  tractions,  which  are  used 

to  close  the  cracks.  Denote  the  solution  as  where  1  indicates  that  this  is  the  solution 

for  the  first  iteration. 

S\EM.  T^'^)  =  + 

2.  Reverse  the  crack  surface  traction  obtained  in  the  previous  step  and  apply  it  as  the  load  on 
the  crack  surfaces  and  solve  the  Pan  a-  Denote  the  solution  as 


3.  Find  the  tractions  on  the  boundary  F,  and  the  displacements  on  the  boundary  F«  from  the 
analytical  solutions  obtained  in  the  previous  step.  Reverse  them  as  the  load  for  Pfem-  Find 
the  crack  closing  tractions  from  the  solution  . 

seem.  = 
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FEAM  converges 
if  the  work  done  in  Pfjj^jTE  is 
less  than  twice  of  the  work  done  in 


the  infinite  body  with  any 
arbitrary  crack  surface 
displacement 

/ 

cxTTTTo:]  r  i 
^i.-L ugfl  \ 

\  \ 

\  IL  \ 


the  finite  body  with  the 
same  crack  surface 
displacement 


(b)  PpiNiTE 


Figure  2.15:  Convergence  Criterion 


4.  Repeat  steps  2  and  3  until  the  residual  load  is  small  enough  to  be  ignored. 


_sANA. 


cFEM. 
•^i+l  • 


for  1  =  2,3,.... 


The  solution  to  the  original  problem  is  the  summation  of  all  those  obtained  in  the  alternating 
procedure,  i.e. 

+  (2.24) 

1=1  '' 


2.5.4  2D  FEAM  for  straight  cracks 

The  analytical  solutions  for  a  crack  in  an  infinite  body,  subjected  to  piecewise  linear  crack  sur¬ 
face  traction  [Wang  and  Atluri(1996)],  is  used  for  the  construction  of  the  2D  FEAM  for  straight 
cracks  in  a  finite  sheet.  Since  the  Finite  Element  Alternating  Method  is  based  on  the  superposition 
principle,  the  complicated  rivet  force  exerted  on  the  skin  by  the  rivets  can  be  easily  taken  into 
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Figure  2.16:  Superposition  Principle  Used  in  the  Finite  Element  Alternating  Method 


account.  Consider  the  local  model  of  the  isolated  cracked  skin  [see  Fig.  2.1 1].  The  skin  is  being 
subjected  to  far-field  tractions,  and  the  stiffener  to  reaction  forces.  The  stress-intensity  factors  for 
single  or  multiple  cracks  (including  Widespread  Fatigue  Damage)  in  the  skin  can  be  determined  in 
the  local  analysis  using  the  Finite  Element  Alternating  Method  (FEAM),  while  still  using  a  coarse 
finite  element  mesh.  The  problem  in  Fig.  2.1 1  can  be  solved  with  the  FEAM  depicted  in  Fig.  2.16, 
wherein  it  can  be  seen  that  the  problem  of  Fig.  2. 1 1 ,  can  be  identified  with  the  problem  labeled  as 
the  “original  problem”  in  Fig.  2.16. 

Essentially,  it  is  a  fixed  point  iteration  scheme  which  solves  the  superposition  of  the  following 
two  problems  [see  Fig.  2.16]. 

1.  the  uncracked,  finite-sized  skin  subjected  to  external  loads  (including  the  reaction  forces 
exerted  by  the  stiffeners  on  the  skin)  and  unknown  external  boundary  loads; 

2.  a  crack  in  an  infinite  sheet  subjected  to  a  crack  surface  traction 

The  crack  surface  traction  in  the  infinite  sheet  cancels  out  the  cohesive  traction  at  the  location  of 
the  crack  in  the  first  problem;  while  the  unknown  external  boundary  loads  in  the  first  problem 
cancel  out  the  tractions  at  the  same  location  in  the  second  problem.  Thus,  the  original  problem, 
i.e.  a  cracked  sheet  of  a  finite  dimension  subjected  to  reaction  forces  exerted  by  the  stiffeners  with 
other  boundary  loadings,  is  exactly  the  superposition  of  these  two  problems. 
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The  FEAM  is  used  to  solve  this  superposition  problem.  First,  reverse  the  cohesive  tractions 
at  the  location  of  the  crack  in  the  finite  element  model  of  the  uncracked  skin  and  use  them  as  the 
load  acting  on  the  crack  in  an  infinite  sheet.  Then,  reverse  the  residuals  at  the  locations  of  the 
far  field  boundries  in  the  infinite  sheet  and  apply  them  as  loads  acting  on  the  boundaries  of  the 
uncracked  skin.  In  this  way,  the  cohesive  tractions  at  the  location  of  the  cracks  and  the  residuals 
at  the  locations  of  the  external  boundaries  are  corrected  iteratively.  This  procedure  converges  very 
fast,  usually  in  two  or  three  iterations.  A  flow  chart  illustrating  the  FEAM  is  shown  in  Fig.  2.17. 

Fracture  mechanics  parameters  can  be  found  accurately  because  the  near  crack  tip  fields  are 
captured  exactly  by  the  analytical  solutions.  Coarser  meshes  can  be  used  in  the  finite  element 
analysis  because  the  cracks  are  not  modeled  explicitly.  The  finite  element  method  is  only  used  to 
compute  the  cohesive  tractions  at  the  crack  location,  which  has  a  smooth  distribution.  Therefore,  a 
very  coarse  mesh  can  be  used.  Fig.  2.18  shows  typical  finite  element  meshes  around  the  crack  tip, 
when  a)  the  Equivalent  Domain  Integral  (EDI)  based  method  is  used  to  evaluate  stress  intensity 
factors;  or,  b)  the  FEAM  is  used.  In  Fig.  2.18,  the  EDI  based  method  also  uses  singular  quarter- 
point  elements.  The  simplicity  of  the  FEAM  mesh  relative  to  the  EDI  mesh,  which  must  explicitly 
model  crack  tips  is  apparent  from  this  figure. 

In  a  parametric  analysis  of  various  crack  sizes,  such  as  is  necessary  in  fatigue  calculations,  the 
stiffness  matrix  of  the  finite  element  model  is  decomposed  only  once,  since  the  stiffness  of  the 
uncracked  structure  remains  the  same  for  all  crack  sizes.  In  the  other  approaches,  such  as  those 
using  singular/hybrid  type  special  crack-tip  finite  elements  or  EDI  methods,  the  cracks  must  be 
modeled  explicitly.  Therefore,  the  global  stiffness  matrix  must  be  computed  and  decomposed  for 
each  crack  size.  Thus,  the  FEAM  is  very  efficient  in  terms  of  both  computational  time  and  human 
effort  (i.e.  mesh  generation)  when  applied  to  problems  such  as  fatigue  crack  growth. 

Finally,  it  is  noted  that  a  simple  superposition  method  can  be  used  to  construct  the  solution 
for  multiple  cracks  in  an  infinite  domain,  subjected  to  arbitrary  crack  surface  tractions,  using  the 
solution  for  a  single  crack  in  an  infinite  domain  (see  the  Appendix  for  theoretical  details).  With 
the  solution  for  multiple  cracks  in  an  infinite  domain,  the  FEAM  can  be  used  to  solve  problems 
of  multiple  cracks  with  arbitrary  crack  lengths  and  orientations  at  arbitrary  locations.  This  is 
particularly  useful  in  the  treatment  of  Widespread  Fatigue  Damage  (WFD). 

2.5.5  3D  FEAM  for  surface  flaws  and  comer  cracks 

The  analytical  solution  for  an  embedded  elliptical  crack  in  an  infinite  body,  subjected  to  arbitray 
crack  face  tractions  [Vijayakumar  and  Atluri  (1981)],  is  used  to  construct  the  3D  FEAM  for  surface 
flaws  and  comer  cracks.  A  20-node  second  order  brick  element  is  used  in  the  3D  FEAM  module. 

The  finite  element  method  is  used  to  analyze  the  uncracked  solid.  Non-zero  stresses  are  calcu¬ 
lated  at  the  location  of  the  actual  crack.  These  stresses  must  be  removed  in  order  to  create  a  traction 
free  crack  as  in  the  actual  problem.  The  infinite  body  with  an  embedded  crack  has  a  solution  which 
is  valid  for  an  arbitrary  distribution  of  tractions  on  the  crack  face.  The  detailed  steps  involved  in 
the  FEAM  for  a  crack  in  a  finite  body  are  as  follows. 

1.  Solve  the  uncracked  finite  body  under  the  given  external  loads  using  the  finite  element 
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(a) 


(b) 


Figure  2.18:  The  Finite  Element  Mesh  When  a)  The  EDI  Method  is  Used;  b)  The  Finite  Element 
Alternating  Method  is  Used 


method.  The  uncracked  body  has  the  same  geometry  as  in  the  given  problem  except  for 
the  crack.  For  example,  when  a  crack  emanates  from  a  hole  in  a  structure,  the  hole  must  still 
be  analyzed  in  the  uncracked  structure. 

2.  Using  the  finite  element  solution,  the  program  computes  the  stresses  at  the  crack  location. 

3.  It  then  compares  the  residual  stresses  calculated  in  Step  2  with  a  permissible  stress  magni¬ 
tude. 

4.  The  residual  stresses  at  the  crack  location  as  computed  in  Step  2  are  reversed  to  create  the 
traction  free  crack  faces  as  in  the  given  problem.  From  this,  the  program  determines  a 
polynomial  form  for  these  stresses  using  a  ’’least  squares  fit”. 

5.  The  analytical  solution  to  the  infinite  body  problem  with  the  crack  subject  to  the  polynomial 
loading  calculated  in  Step  4  is  now  obtained. 
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6.  The  stress  intensity  factors  for  the  current  iteration  are  then  calculated  from  the  analytical 
solution. 

7.  The  residual  stresses  on  the  external  surfaces  of  the  body  due  to  the  applied  loads  on  the 
crack  faces,  are  now  computed.  To  satisfy  the  given  traction  boundary  conditions  at  the 
external  boundaries,  the  residual  stresses  on  the  external  surfaces  of  the  body  are  reversed 
and  this  allows  calculation  of  the  equivalent  nodal  forces. 

8.  Consider  the  nodal  forces  in  Step  7  as  externally  applied  loads  acting  on  the  uncracked  body. 


All  the  steps  in  the  iteration  process  are  repeated  until  the  residual  stresses  on  the  crack  surface 
become  negligible.  It  has  been  observed  that  this  iteration  process  typically  takes  three  or  four 
steps.  The  overall  stress  intensity  factor  solution  is  obtained  by  adding  the  stress  intensity  factor 
solutions  for  all  iterations. 

A  recent  development  in  the  3D  FEAM  during  Phase  II  of  this  project  is  the  development  of  a 
numerical  scheme  for  the  evaluation  of  the  analytical  solution  for  circular  or  near  circular  cracks. 
The  numerical  procedure  developed  by  Nishioka  and  Atluri(1983)  becomes  numerically  unstable 
as  the  aspect  ratio  of  the  ellipse  tends  to  1.  When  the  ellipse  becomes  a  circle,  this  numerical 
procedure  is  no  longer  valid.  However,  the  analytical  solution  based  on  the  ellipsoidal  potential 
[Vijayakumar  and  Atluri  (1981)]  is  still  valid.  An  alternative  numerical  procedure  was  developed 
during  Phase  II  of  this  project.  This  alternative  numerical  procedure  is  documented  in  this  section; 
while  the  reader  is  refered  to  Vijayakumar  and  Atluri  (1981)  and  Nishioka  and  Atluri(1983)  for 
other  details  on  the  analytical  solution  based  on  the  ellipsoidal  potentials. 

Numerical  difficulties  arise  in  the  evalution  of  the  generic  elliptical  integral  Lijm,  defined  as 


■^ijm 


=/“ 
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sn  nd  nc  du 


(2.25) 


where  nd  u,  sn  u,  and  nc  u  are  Jacobian  elliptical  functions. 

The  numerical  procedure  introduced  by  Nishioka  and  Atluri(1983)  uses  a  recursive  formula 
starting  from  the  evalution  of  the  integrals 


r  - 
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n  =  0, ihl, it2, . . . 


(2.26) 


Let  a  be  the  semi-major  axis;  b  be  the  semi-minor  axis;  and  k'  =^bl a  be  the  aspect  ratio  of 
the  ellipse.  =  1  —  k'^.  It  is  seen  that  as  k'  1,  >  0.  The  recursive  formula  makes  use  of  the 

identity  k^sn ^u  +  dn^u=  I  to  evaluate  sn  for  a  given  dn ^u.  However,  this  formula  breaks  down 
when  ^  =  0,  since  dn^M  =  1  in  this  case.  Alternatively,  we  can  evaluate  dn^M  for  a  given  sn^M. 
Thus,  we  can  start  the  recursive  formula  from  the  evaluation  of  the  integrals 

ru 

/;„=  /  sn^'^udu  «  =  0,1,2,...,  (2.27) 

u  0 

Two  numerical  schemes  were  developed  to  evaluate  in  Phase  II  of  this  project.  One  is  for 
the  case  where  k  is  close  to  1 ;  and  the  other  is  used  for  small  k. 
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For  small  k,  the  integral  is  represented  in  terms  of  a  series. 

4"„=  [  sn%</M=  n  =  0,l,2,...,M+2 
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where  sine])  =  sn  u,  and 
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It  is  evaluated  using  the  following  recursive  formula. 


where  =  <|)  and 
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For  large  k,  the  recursive  formula  is 


jn+\ _  ^n^sn 

CM 

Pn 


^sn 


where 


Pn  =  {l  +  2n)k^,  qn  =  2n{l  +  k^),  r„  =  (l-2n). 
a„  =  a„_i  X  sn^M,  ai  =  sn  «  cn  m  dn  u 

ro-FU)  ji  m-m 


(2.28) 


(2.29) 


(2.30) 


(2.31) 


(2.32) 

(2.33) 

(2.34) 

(2.35) 


Once  /"  are  evaluated,  Lum  can  be  computed  using  the  following  scheme.  First,  compute 
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<li  =  2[(2;  +  m  -  /)  +  {i  -  j)k^\ 

(2.42) 

0  =  l  +  2(/-y-m) 

(2.43) 

tty  =  ay_i  X  nd^M 

(2.44) 

tto  =  — u  cn  '“^'"w/dn  u 

(2.45) 

Here,  m  =  Oorm=— 1.  Therefore,  no  singular  term  exists. 

Then,  compute  Lij^  for  m  =  0, 1 , . . (M—  1  —  j)  using  the  following  recursive  formula. 

^mLi  jm  +  ^ pnA‘,7,(m- 1)  + 
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where 
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(2.49) 
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(2.51) 

Singular  terms  exist  as  cn  m  ->•  0.  In  this  case,  sn  w  =  1  and  dn  =  Vl  - 

r  2m+lr  sn^'+^Mnd^-^~^ 

lim  cn^'"+'L,.  = - 

cntt->0  hJ,(>n+i) 

(2.52) 

The  recursive  formula  can  be  derived  from  the  integration  of 

[sn  *  M  nd  *  M  nc  ^  m]  ^ 

Note,  the  following  recursive  formula  is  only  used  for  the  case  where  j 

=  0,  m  =  0  and  k  is 

large  enough. 
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where 

Pi=  [l+2{i-j-m)]k^ 

(2.54) 

qi  =  2[{i-m)  +  (i-j)k'^] 

(2.55) 
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(2.56) 
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(2.58) 

One  recent  application  [O'Donoghue,  Atluri  and  Pipkins  (1995)]  of  the  FEAM  was  to  the 
analysis  of  fatigue  cracking  in  the  lower  wing  skin  of  the  U.S.  Air  Force  C-141B  (Fig.  2.19).  In 
this  application,  the  growth  of  corner  cracks  from  the  weep  holes  located  in  the  integral  risers  of 
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Figure  2.19:  U.S.  Air  Force  C-141B 


2.20:  Cut-Out  Lower  Wing  Panel  from  the  C-141B  Showing  Weep 


Figure  2.21 :  Cross-Section  of  Failed  Lower  Wing  Panel  of  the  C- 14  IB 


the  wing  skin  was  modeled  with  the  FEAM  (Fig.  2.20  and  Fig.  2.21).  Fatigue  crack  growth 
predictions  were  made  wherein  the  FEAM  was  used  to  generate  stress  intensity  factors  which 
were  intum  used  in  a  Forman  equation.  The  load  spectrum  used  was  comprised  of  peak-valley 
pairs  representing  3027  equivalent  flight  hours.  Comparisons  with  limited  test  data  showed  good 
correlation  between  the  FEAM  based  fatigue  crack  growth  predictions  and  the  experimental  data. 

2.5.6  Distributed-Dislocation-base  FEAM  for  curved  cracks 

A  curved  crack  can  be  treated  as  a  distribution  of  dislocations.  To  construct  the  analytical  solution 
for  a  curved  crack  in  an  infinite  sheet,  one  can  solve  for  the  dislocation  density  for  the  curved 
crack.  Using  a  complex  stress  function  approach,  the  dislocation  density  can  be  related  to  the 
crack  surface  traction  [Park  and  Atluri(1998),  Chen(1993),  and  Cheung  and  Chen(1987)].  Once 
the  dislocation  density  is  solved  for,  the  stress  and  displacement  in  the  infinite  sheet,  subjected  to 
the  prescribed  crack  surface  traction,  can  be  obtained. 

The  formulation  in  Park  and  Atluri(1998)  is  used  for  the  construction  of  the  dislocation-based 
analytical  solutions  for  a  curved  crack,  subjected  to  arbitrary  crack  surface  tractions.  Using  these 
solutions  with  the  finite  element  solution  for  a  sheet  of  finite  size,  a  Dislocation-based  FEAM  is 
implemented  for  the  curved  cracks  in  a  finite  sheet. 

2.6  Fatigue  crack  growth 


Fatigue  Crack  Growth 

The  problem  of  fatigue  crack  growth  is  of  considerable  practical  importance  when  designing  a 
structure  which  satisfies  the  DTR.  To  successfully  employ  damage  tolerance  principles,  an  accurate 
determination  is  required  of  the  number  of  load  cycles  to  failure  in  a  component.  These  estimates 
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will  have  a  significant  influence  in  the  design  and  maintenance  of  a  safe  structure  such  as  in  the 
scheduling  of  inspection  intervals. 

Fatigue  crack  growth  frequently  occurs  when  a  flawed  component  is  subject  to  some  form  of 
cyclic  loading.  Here  the  crack  growth  is  termed  subcritical  since,  due  to  the  cyclic  loading,  it  takes 
place  at  stress  intensity  factor  levels  that  are  less  than  the  fracture  toughness  of  the  material.  The 
form  of  the  cyclic  loading  is  also  of  great  importance  such  as  whether  it  has  constant  amplitude  or 
has  a  variable  amplitude.  The  damage  tolerance  module  has  the  capability  to  model  fatigue  crack 
under  conditions  of  constant  amplitude  and  variable  amplitude  loading.  This  capability  is  limited 
to  self  similar  (i.e.  Mode  I)  crack  growth. 

The  crack  growth  calculations  will  be  performed  based  on  stress  intensity  factors  obtained 
either  by  the  FEAM  or  user  supplied  beta  factors.  For  reasons  of  computational  efficiency  and 
accuracy,  the  user  can  use  the  FEAM  to  update  the  beta  factor  after  the  crack  growth  exceeds  a 
specific  amount.  Fatigue  crack  growth  is  carried  out  using  the  beta  factor,  assuming  that  a  small 
amount  of  crack  growth  does  not  the  effect  the  beta  factor.  When  the  amount  of  crack  growth 
exceeds  the  threshold  specified  by  the  user,  the  FEAM  can  be  invoked  to  update  the  beta  factor. 

The  FEAM  has  made  the  use  of  finite  elements  in  fatigue  calculations  feasible.  The  reason 
being  that  the  global  stiffness  matrix  of  the  finite  element  model  is  assembled  and  decomposed  only 
once,  since  the  stiffness  of  the  uncracked  structure  remains  the  same  for  all  crack  sizes.  In  other 
finite  element  approaches,  such  as  those  using  singular/hybrid  type  special  crack-tip  finite  elements 
or  using  EDI  methods,  the  cracks  must  be  modeled  explicitly.  Therefore,  the  global  stiffness  matrix 
must  be  assembled  and  decomposed  for  each  crack  size.  The  assembly  and  decomposition  of  the 
global  stiffness  matrix  accounts  for  about  80%  of  the  computational  time  required  in  a  typical  finite 
element  analysis.  Given  that  the  FEAM  has  to  perform  this  operation  only  once  during  a  fatigue 
analysis,  the  benefits  of  this  approach  over  other  finite  element  techniques  for  fatigue  crack  growth 
are  readily  apparent. 

Load  spectra,  in  terms  of  ASTROS'  load  cases,  is  provided  to  the  damage  tolerance  module 
by  the  USAGE  module  described  in  Nees  (1995).  Since  the  USAGE  module  defines  aircraft  ma¬ 
neuvers  in  terms  of  ASTROS  load  cases,  the  load  spectra  at  any  point  in  the  aircraft  structure  is 
automatically  determined  by  extracting  stresses  at  that  point  from  the  ASTROS  solution  database. 
Use  is  made  of  USAGE  features  such  as  repeating  and  blocking  of  data  to  minimize  storage  and/or 
computational  requirements  of  variable  amplitude  crack  growth  calculations. 

Numerous  studies  have  been  conducted  on  the  characteristics  of  fatigue  crack  growth.  It  has 
been  established  that  when  the  plastic  or  inelastic  zone  in  the  vicinity  of  the  crack  is  small,  then  the 
stress  intensity  factor  is  the  governing  parameter  during  crack  growth  [Paris  and  Erdogan  (1963)]. 
In  general,  the  crack  growth  rate  is  a  function  of  stress  intensity  factor  change,  which  is  given  by; 

AK  =  Kmax-Krr,in  (2.59) 

where  K^ax  is  the  maximum  stress  intensity  factor  during  the  load  cycle  and  K^m  is  the  minimum 
value  of  the  stress  intensity  factor.  Based  on  the  minimum  and  maximum  loads,  it  is  customary  to 
define  the  parameter  R  where: 

R^KminlKnu^  (2-60) 
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An  important  point  must  be  made  in  relation  to  the  functional  relationship  between  the  fatigue 
crack  growth  rate  and  the  change  in  stress  intensity  factor.  This  crack  growth  rate  function  can  be 
partitioned  into  three  separate  regions.  At  low  values  of  6K,  there  is  very  little  crack  growth  with  a 
negligible  crack  growth  rate.  Therefore,  it  can  be  stated  that  there  exists  a  A/sT  below  which  there  is 
no  crack  growth.  This  quantity  is  referred  to  as  the  threshold  stress  intensity  factor  and  is  denoted 
as  [AK)(fi.  At  higher  values  of  AK,  crack  growth  takes  place.  It  has  been  observed  experimentally 
that  this  curve,  relating  the  crack  growth  rate,  da/dN,  to  AK,  is  usually  linear  on  a  log-log  plot  and 
this  corresponds  to  a  power  law  relation  between  daj dN  and  AK.  This  is  commonly  referred  to  as 
the  Paris  relation  and  is  given  as  [Paris,  Gomez  and  Anderson  (1961)]: 

^  =  C(A«)"  (2.61) 

where  a  is  the  crack  length  and  N  is  the  number  of  load  cycles.  The  quantities  C  and  n  are  material 
dependent  constants.  At  higher  values  of  AK  the  stress  intensity  factor  is  approaching  the  fracture 
toughness  of  the  material,  Kc.  The  crack  growth  rate  will  increase  significantly,  eventually  leading 
to  the  onset  of  rapid  unstable  crack  growth. 

Recognizing  that  several  distinct  phases  in  fatigue  crack  growth  exist,  a  more  general  form  for 
the  relationship  between  crack  growth  rate  and  stress  intensity  factor  is  expressed  as: 

dN  { ( 1  -  R)Kc  -  AKY 

where  m,  p  and  q  are  constants  that  relate  to  the  particular  crack  growth  relation  that  is  being  used. 
By  assigning  different  values  to  these  quantities  some  of  the  well  known  crack  growth  relations 
can  be  recovered.  For  example,  when  m  =  p  =  ^  =  0,  the  Paris  relation  is  obtained.  The  Forman 
relation  [Forman,  Kearney  and  Engle  (1967)],  which  accounts  for  high  crack  growth  rates  and 
instability,  is  recovered  by  setting  m  =  p  =  0  and  ^  =  1.  By  setting  =  {M^  -l)n, 

the  Walker  relation  [Walker  (1970)]  is  derived,  where  is  an  exponent  in  the  Walker  relation. 

With  stress  intensity  factor  solutions  and  a  crack  growth  relation  in  hand,  the  ultimate  objective 
of  fatigue  crack  growth  calculations,  to  calculate  the  number  of  cycles  for  a  crack  to  grow  by  a 
specified  amount,  can  be  carried  out.  This  is  done  by  integrating  the  growth  relation  [Equation 
(2.61)  or  (2.62)].  For  example,  the  equation  (2.61)  is  integrated  as  follows.  Sinee 

AK  =  Acy/na 


Eq.  2.61  can  be  rewritten  as 


n/2 


where  Ci  —  C  [PAa/tt]  .  Assuming  that  the  beta  factor  does  not  ehange  for  the  small  amount  of 
crack  growth  in  N  cycles,  the  amount  of  crack  growth  is 
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or 


P[aa:)^  1 

""  ^  —Uo  for  n  =  2 

where  ao  is  the  half  crack  length  before  the  N  cycles  crack  growth.  The  integration  is  carried  out 
for  a  specific  design  life  to  find  the  final  crack  length  at  the  end  of  the  end  of  the  design  life  of  the 
aircraft.  All  of  these  steps  will  be  carried  out  automatically  by  the  damage  tolerance  module. 

Self-similar  growth  of  through  skin  cracks  is  straightforward  to  carry  out  using  the  above  pro¬ 
cedure.  However,  for  an  elliptical  or  part  elliptical  crack  since  the  stress  intensity  factor  distribution 
is  generally  not  constant  along  the  elliptical  crack  front,  the  crack  growth  rates  will  not  be  the  same 
in  every  direction.  Consequently,  the  shape  of  the  crack  will  change.  It  has  been  often  observed  in 
practice  that  a  semi-elliptical  surface  crack,  initially  having  a  small  aspect  ratio,  will  have  a  larger 
crack  growth  rate  in  the  minor  axis  direction.  The  damage  tolerance  module  computes  the  amount 
of  crack  growth  for  the  major  and  minor  axis,  assuming  that  the  orientation  of  the  major/minor 
axis  does  not  change.  The  situation  where  the  minor  axis  grows  faster  than  the  major  axis  so  that 
it  becomes  the  major  axis  after  a  numbers  of  cycles  is  allowed  by  the  damage  tolerance  module. 

2.7  Optimization 

The  principal  objective  of  multi-disciplinary  structural  optimization  is  to  minimize  an  objective 
function  such  as  structural  weight  and/or  cost  subject  to  discipline  specific  constraints  (strength, 
flutter,  stiffness,  etc.).  ASTROS  provides  a  sophisticated  software  framework  for  performing  such 
optimizations.  ASTROS  supports  both  preliminary  design  and  design  modifications  that  occur 
later  in  the  product  life  cycle.  It  combines  finite  element  modeling  and  analysis  techniques  with 
efficient  optimization  algorithms  to  deliver  significant  reductions  in  the  time  required  to  develop 
superior  designs  of  structures.  The  finite  element  modeling  is  based  on  the  standard  NASTRAN 
bulk  data  format. 

The  preliminary  design  stage  is  where  ASTROS  capabilities  can  be  used  to  the  fullest.  Typ¬ 
ically,  the  aerodynamic  configuration,  materials  and  design  conditions  have  been  defined  at  this 
stage.  The  preliminary  design  challenge  is  then  to  determine  the  optimal  (i.e.  least  weight  and/or 
cost)  structural  configuration  which  satisfies  constraints  imposed  by  multiple  engineering  disci¬ 
plines.  ASTROS  supports  a  wide  variety  of  constraints.  These  include:  Tsai-Wu  stress  crite¬ 
ria,  von-Mises  stress;  stiffness;  natural  frequency;  flutter;  laminate  composition;  panel  and  beam 
buckling;  and  aeroelastic  lift  and  control  effectiveness.  With  the  addition  of  the  damage  tolerance 
module,  residual  strength  and  residual  life  based  constraints  are  now  part  of  this  list.  ASTROS 
performs  the  optimization  task  by  automatically  changing  design  variables.  The  design  variables 
which  ASTROS  uses  are  related  to  properties  of  the  finite  elements  used  to  model  the  structure 
being  optimized.  The  design  variables  are  the  thickness  of  shell  elements  and  cross-sectional  area 
of  beam  elements.  For  laminates,  the  thickness  of  individual  plys  is  the  design  variable.  In  order  to 
insure  that  the  optimized  structure  is  realistic  from  a  manufacturing  point  of  view  ASTROS  sup¬ 
ports  design  variable  linking.  Design  variable  linking  allows  the  grouping  of  elements  for  constant 
thickness  structure  or  the  definition  of  shape  functions  for  tapered  (thickness)  structure. 
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The  damage  tolerance  module  evaluates  the  residual  strength  and  residual  life  based  constraints 
based  on  damage  tolerance  analysis.  However,  automated  redesign  based  on  the  damage  tolerance 
constraint  is  still  not  well  resolved.  Based  on  the  factor  that  damage  tolerance  analysis  takes  a 
significant  amount  of  time,  while  the  damage  tolerance  constraint  is  not  the  most  important  con¬ 
straint  during  the  preliminary  design,  it  is  suggested  that  the  user  used  beta-factor  based  automated 
redesign  strategy.  The  details  follow.  The  user  can  use  the  damage  tolerance  module  to  evaluate 
the  beta  factors  as  correction  factors  for  an  analytical  solution  for  a  simple  case.  The  constraint 
can  than  be  expressed  explicitly  in  using  a  formula  involving  the  beta  factor.  This  constraint  can 
be  defined  by  the  user  in  terms  of  a  user  defined  function  in  ASTROS.  Assuming  that  the  beta 
factor  does  not  change  during  the  redesign,  one  can  use  the  existing  capability  of  ASTROS  to  per¬ 
form  automated  redesign.  After  the  redesign,  the  beta  factor  can  be  re-evaluated  to  ensure  that  the 
damage  tolerance  constraint  is  met.  This  approach  can  significantly  reduce  the  computational  time 
associated  with  the  computation  of  damage  tolerance  constraints. 

2.8  Interfacing  with  ASTROS 

This  section  briefly  summaries  the  interfacing  techniques  used  to  integrate  the  damage  tolerance 
module  with  ASTROS  (Automated  STRuctural  Optimization  System).  More  details  can  be  found 
in  the  Interface  Design  Document  for  this  project. 

ASTROS  is  an  extendible  system  built  on  top  of  an  engineering  database  and  the  Matrix 
Analysis  Problem  Oriented  Language  (MAPOL).  The  database  provides  a  channel  for  the  inter¬ 
module  communication  in  ASTROS;  while  MAPOL  provides  extendibility  to  ASTROS.  During 
the  integration  of  the  damage  tolerance  module  into  ASTROS,  an  extremely  powerful  “gluing  tool” 
-  Tool  Command  Language  (TCL)  is  introduced  into  ASTROS. 

The  damage  tolerance  module  is  a  complex  system  by  itself.  It  includes  customized  mesh  gen¬ 
erators,  finite  element  alternating  [Schwartz-Neumann  Alternating]  codes  (for  2D  straight  cracks, 
2D  curved  cracks  and  3D  elliptical/circular  cracks),  an  automated  Global/Local  analyzer,  and  a 
buckling  analyzer.  They  were  developed  using  different  computer  languages,  such  as  C,  C-H-  and 
FORTRAN.  In  order  to  integrate  the  damage  tolerance  module  into  ASTROS,  a  TCL  interpreter  is 
planted  into  ASTROS  as  a  module  invokable  from  a  MAPOL  sequence. 

Fig.  2.22  shows  the  interfacing  strategy  for  incorperating  damage  tolerance  into  ASTROS. 
It  is  seen  in  Fig.  2.22  that  the  Damage  Tolerance  Module  communicates  with  the  other  Built-in 
Modules  of  ASTROS  through  the  Engineering  Database.  The  extended  modules  are  controlled  by 
TCL  scripts,  which  are  interpreted  by  the  embedded  TCL  interpreter.  The  interpreter  is  in  turn 
controlled  by  the  modified  MAPOL  sequence. 

Interfacing  aspects,  related  to  MAPOL,  DATABASE  and  TCL,  are  documented  in  this  chapter. 
Database  entities  used  by  the  damage  tolerance  module  to  communicate  with  other  modules  are 
described  in  the  Interface  Design  Document. 
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2.8.1  MAPOL 


As  seen  in  Fig.  2.22,  the  MAPOL  sequence  does  not  control  the  damage  tolerance  module  directly. 

It  interacts  with  the  TCL  module  only. 

A  single  argument  (INTEGER)  is  passed  to  the  TCL  module  from  the  MAPOL  calling  sequence. 
This  argument  indicates  the  position  in  the  MAPOL  sequence,  from  which  the  MAPOL  call  is 
placed.  This  status  indicator  is  set  as  the  value  of  the  global  TCL  variable  astrosStatelD.  De¬ 
pending  on  the  value  of  astrosStatelD,  the  TCL  module  i)  creates  an  interpreter;  ii)  performs 
input  data  check;  iii)  performs  damage  tolerance  analysis;  or  iv)  deletes  the  interpreter,  etc.  No 
other  direct  data  communication  is  made  between  the  MAPOL  calling  sequence  and  the  damage 
tolerance  module.  The  damage  tolerance  module  communicates  with  the  rest  of  the  modules  in 
ASTROS,  though  the  database  entities,  to  get  the  data  required  for  damage  tolerance  analyses.  It 
also  stores  the  results  in  the  database  for  other  modules  to  access. 

2.8.2  DATABASE 

ASTROS  relies  on  the  engineering  database  for  inter-module  conununication.  Originally  based  on 
the  Computer  Automated  Design  Database  (CADDB),  which  was  the  heart  of  ASTROS  during  its 
development,  the  commercially  supported  ASTROS  is  now  based  on  the  eBASE  database  devel¬ 
oped  by  UAL  In  this  project,  the  set  of  interface  routines  for  CADDB  are  used  to  access  database 
entries.  To  facilitate  the  migration  from  CADDB  to  eBASE,  the  damage  tolerance  module  uses 
a  set  of  TCL  commands  to  access  the  database  via  the  CADDB  interface.  To  migrate  from  the 
CADDB  interface  to  the  eBASE  interface,  it  will  be  necessary  to  re-implement  only  those  TCL 
commands  using  the  eBASE  access  routines. 

2.8.3  TCL 

Tool  Command  Language  (TCL)  is  a  simple,  yet  robust  scripting  language.  It  is  an  excellent  script¬ 
ing  language  for  extending  the  functionality  of  existing  programs.  TCL  was  originally  developed 
at  the  University  of  California,  Berkeley.  The  development  was  later  shifted  to  Sun  Microsystems, 

Inc.  Currently,  it  is  commercially  supported  by  Scriptics,  Inc.  The  core  system  of  TCL  is  freely 
distributed  with  the  source  code.  Although  it  is  distributed  without  a  charge,  it  is  very  stable  and  of 
commercial  quality.  TCL  has  become  an  increasingly  popular  cross-platform  scripting  language. 
Many  commercial  software  have  been  developed  based  on  TCL. 

TCL  consists  of  a  scripting  language  and  an  interpreter  for  that  language.  TCL  interpreter  is 
designed  such  that  it  can  be  easily  embedded  into  other  applications.  As  a  language,  it  is  much  like 
the  UNIX  shell  language.  There  is  very  little  syntax;  and  it  is  very  easy  to  learn.  TCL  has  been 
used  to  assemble  software  modules  that  have  been  built  in  system  programming  languages  like  C, 
C-H-  and  FORTRAN.  These  building  blocks  appear  as  commands,  or  verbs,  in  TCL. 

Online  user  manuals  are  available  at  Scriptics'  web  site  (http ;  / / www .  scriptics .  com/man 
/tcl8 . 0/contents .  htm).  Much  more  helpful  information  can  be  found  at  http :  /  /  www .  scriptics . 
com/resource/. 
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The  damage  tolerance  module  is  actually  a  complex  system.  It  includes  customized  mesh  gen¬ 
erators,  finite  element  alternating  codes  [based  on  the  Schwartz-Neumann  alternating  method]  (for 
2D  straight  cracks,  2D  curved  cracks  and  3D  elliptical/circular  cracks),  an  automated  Global/Local 
analyzer,  and  a  buckling  analyzer.  They  were  developed  using  different  computer  languages,  such 
as  C,  C++  and  FORTRAN.  In  order  to  integrate  the  damage  tolerance  module  into  ASTROS,  a 
TCL  interpreter  was  planted  into  ASTROS  as  a  module  invokable  from  the  MAPOL  sequence. 
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Engineering  Database:  CADDB 


Figure  2.22:  Integration  of  damage  tolerance  module  with  ASTROS 
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CHAPTER  III 

SUMMARY  OF  PHASE  II  ACCOMPLISHMENTS 


A  summary  of  the  Phase  II  accomplishments  in  this  SBIR  project  is  presented  in  this  chapter. 

During  Phase  II,  a  damage  tolerance  module  is  implemented  and  integrated  with  ASTROS. 
To  integrate  the  damage  tolerance  module  into  ASTROS,  the  TCL  (Tool  Command  Language) 
interpreter  has  been  implanted  into  ASTROS  as  a  gluing  tool. 

Customized  damage  tolerance  models  have  been  implemented  so  that  the  user  can  easily  per¬ 
form  these  types  of  damage  tolerance  analysis  with  very  little  effort.  A  master  element  approach 
has  been  used  to  minimize  the  impact  of  damage  tolerance  analysis  on  the  data  preparation  for  the 
preliminary  design  model.  With  a  few  extra  bulk  data  cards,  the  user  can  perform  damage  tolerance 
analysis. 

The  customized  damage  tolerance  models  that  have  been  implemented  are: 

•  Discrete  Source  Damage  Model:  a  lead  crack  in  a  stiffened  panel  with/without  the  presence 
of  a  broken  central  stiffener  [Fig.  2.7] 

•  BuckDel  Model:  buckling  of  a  composite  panel  in  the  presence  of  a  delamination  [Fig.  2.8] 

•  Straight  Crack  Model:  a  panel  with  a  centered  crack  [Fig.  2.1] 

•  Rivet  Hole  Crack  model:  one  (or  two)  crack(s)  emanating  from  one  side  (or  both  sides)  of  a 
rivet  hole  [Fig.  2.2] 

•  Curved  Crack  model:  a  panel  with  a  curved  crack  [Fig.  2.5] 

•  Rivet  Hole  Curved  Crack  model:  one  (or  two)  curved  crack(s)  emanating  from  one  side  (or 
both  sides)  of  a  rivet  hole  [Fig.  2.6] 

•  Surface  crack  model:  one  centered  surface  crack  in  a  plate  [Fig.  2.3] 

•  Rivet  Hole  Corner  Crack  model:  two  corner  cracks  emanating  from  both  side  of  a  straight- 
shank  rivet  hole.  [Fig.  2.4] 

An  automated  global-local  analyzer  was  developed  for  the  analysis  of  discrete  source  damage. 
The  automated  global-local  analyzer,  based  on  the  feature  modeling  technique,  has  an  integrated 
geometry  modeler  and  mesh  generator.  It  completely  automates  the  model  simplification,  model 
generation  and  submodel  creation  for  a  hierarchical  analysis. 

A  2D  FEAM  (Finite  element  alternating  method),  3D  FEAM,  and  distributed  dislocation  based 
FEAM  have  been  implemented  in  the  damage  tolerance  module  for  the  analyses  of  2D  straight 
cracks,  3D  surface  flaws/corner  crack,  and  2D  curved  cracks.  An  alternative  numerical  scheme  has 
been  developed  to  overcome  the  difficulties  associated  with  the  cracks  of  a  circular  shape. 
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The  damage  tolerance  module  can  generate  a  cyclic  loading  spectrum,  based  on  the  user  defined 
loading  spectrum  (in  terms  of  ASTROS  loading  cases).  Using  this  cyclic  loading  spectrum,  the 
damage  tolerance  module  can  perform  self-similar  (Mode  I)  fatigue  crack  growth  analysis. 

The  damage  tolerance  module  was  integrated  into  ASTROS  v20p0  on  AIX  4.2.  It  has  been 
ported  on  the  IRIX  operating  system. 
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